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ABSTRACT 


A laboratory spacecraft simulator testbed is first introduced to examine the 
problem of multiple spacecraft interacting in close proximity. This testbed enables 
validation of guidance, navigation and control (GNC) algorithms by combining 6- 
Degrees of Freedom (DoF) computer simulation with 3-DoF Hardware-In-the-Loop 
(HIL) experimentation. The presented 3-DoF spacecraft simulator employs a novel 
control actuator configuration consisting of a Miniature Single Gimbaled Control 
Moment Gyroscope (MSGCMG) and dual on/off cold-gas in-plane vectorable thrusters. 
The dual vectorable thruster design enables simultaneous translation and attitude control 
allowing it to act both in conjunction with the MSGCMG as well as provide sole actuator 
control throughout a commanded closed-path maneuver. Small-time local controllability 
(STLC) of this uniquely actuated system via Lie Algebra methods is formally 
demonstrated and results of experiments conducted on the described testbed are included. 
From this study in 3-DoF, a 6-DoF minimally control actuated asymmetric spacecraft 
design is proposed. Six-DoF control of this underactuated mechanical system is achieved 
via two oppositely mounted hemispherically vectorable thrusters. In order to capitalize 
on the unique nature of this system with only two control torques, a quaternion feedback 
regulator is developed to yield three-axis stabilization of its attitude. This regulator 
capitalizes on recent advancements in generalized inversion and perturbed feedback 
linearizing control to stabilize the dynamics of an underactuated asymmetric spacecraft 
and extends this to include stabilization of the kinematics of the system. Two control 
design methodologies are derived. The first is Lyapunov based, yielding a globally stable 
system, while the second yields local stability within a domain of attraction through 
perturbed feedback linearization. Results of several numerical simulations are presented 
for an asymmetric spacecraft with two bounded body-fixed control torques. The 
proposed attitude control method is not intended to provide attitude maintenance for 
attitude tracking or in the presence of relatively large disturbance torques; however, it 
may prove widely applicable for detumbling and reorientation maneuvers of spacecraft 
with only two available control torques. 
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I. 


INTRODUCTION 


The traditional spacecraft system is a monolithic structure with a single mission 
focused design and lengthy production and qualification schedules coupled with 
enormous cost. Additionally, there rarely, if ever, are any designed preventive 
maintenance plans or re-fueling capability. There has been much research in recent years 
into alternative options. One option involves autonomous on-orbit servicing of current or 
future monolithic spacecraft systems. The U.S. Department of Defense (DoD) embarked 
on a highly successful venture to prove such a concept with the Defense Advanced 
Research Projects Agency’s (DARPA) Orbital Express program. Orbital Express 
demonstrated all of the enabling technologies required for autonomous on-orbit servicing 
to include refueling, component transfer, autonomous satellite grappling and berthing, 
rendezvous, inspection, proximity operations, docking and undocking, and autonomous 
fault recognition and anomaly handling (Kennedy 2008). Another potential option 
involves a paradigm shift from the monolithic spacecraft system to one involving 
multiple interacting spacecraft that can autonomously assemble and reconfigure. 
Numerous benefits are associated with autonomous spacecraft assemblies, ranging from a 
removal of significant intra-modular reliance that provides for parallel design, 
fabrication, assembly and validation processes to the inherent smaller nature of 
fractionated systems that allows for each module to be placed into orbit separately on 
more affordable launch platforms (Mathieu and Weigel 2005). 

A. GROUND-BASED HARDWARE-IN-THE-LOOP EXPERIMENTAL 

VALIDATION 

With respect specifically to the validation process, the significantly reduced 
dimensions and mass of aggregated spacecraft when compared to the traditional 
monolithic spacecraft allow for not only component but even full-scale on-the-ground 
Hardware-In-the-Loop (HIL) experimentation. Likewise, much of the HIL 
experimentation required for on-orbit servicing of traditional spacecraft systems can also 
be accomplished in ground-based laboratories (Creamer 2007). This type of HIL 
experimentation complements analytical methods and numerical simulations by 
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providing a low-risk, relatively low-cost and potentially high-return method for 
validating the technology, navigation techniques and control approaches associated with 
spacecraft systems. Several approaches exist for the actual HIL testing in a laboratory 
environment with respect to spacecraft guidance, navigation and control. One such 
method involves reproduction of the kinematics and vehicle dynamics for 3-DoF (two 
horizontal translational degrees and one rotational degree about the vertical axis) through 
the use of robotic spacecraft simulators that float via planar air bearings on a flat 
horizontal floor. This particular method is currently being employed by several research 
institutions and is the validation method of choice for our research into GNC algorithms 
for proximity operations at the Naval Postgraduate School (Machida, Toda, and Iwata 
1992; Ullman 1993; Corrazzini and How 1998; Marchesi, Angrilli and Venezia 2000; 
Ledebuhr et al. 2001; Nolet, Kong, and Miller 2005; LeMaster, Schaechter, and 
Carrington 2006; Romano, Friedman, and Shay 2007). 

For spacecraft involved in proximity operations, the in-plane and cross-track 
dynamics are decoupled, as modeled by the Hill-Clohessy-Wiltshire (HCW) equations, 
thus the reduction to 3-Degree of Freedom (DoF) does not appear to be a critical limiter. 
One consideration involves the reduction of the vehicle dynamics to one of a double 
integrator. However, the orbital dynamics can be considered to be a disturbance that 
needs to be compensated for by the spacecraft navigation and control system during the 
proximity navigation and assembly phase of multiple systems. Thus, the flat floor 
testbed can be used to capture many of the critical aspects of an actual autonomous 
proximity maneuver that can then be used for validation of numerical simulations. 
Portions of the here-in described testbed, combined with the first generation robotic 
spacecraft simulator of the Spacecraft Robotics Laboratory (SRL) at the Naval 
Postgraduate School (NPS), have been employed to propose and experimentally validate 
control algorithms. The interested reader is referred to Romano et al. (2007) for a full 
description of this robotic spacecraft simulator and the associated HIL experiments 
involving its demonstration of successful autonomous spacecraft approach and docking 
maneuvers to a collaborative target with a prototype docking interface of the Orbital 
Express program. 
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A unique control problem exists, given a requirement for spacecraft aggregates to 
rendezvous and dock during the final phases of assembly and a desire to maximize the 
useable surface area of the spacecraft for power generation, sensor packages, docking 
mechanisms and payloads while minimizing thruster impingement. In fact, control of 
such systems using the standard control actuator configuration of fixed thrusters on each 
face coupled with momentum exchange devices can be challenging, if not impossible. 
For such systems, a new and unique configuration is proposed consisting of vectorable 
thrusters that may capitalize, for instance, on the recently developed carpal robotic joint 
invented by Canfield and Reinholtz (1998) with its hemispherical vector space. It is here 
demonstrated through Lie algebra analytical methods and experimental results that two 
vectorable in-plane thrusters in an opposing configuration can yield a minimum set of 
actuators for a controllable reduced order (3-DoF) system. It will be shown that by 
coupling the proposed set of vectorable in-plane thrusters with a single degree of freedom 
torquer such as a Control Moment Gyroscope, an additional degree of redundancy can be 
gained for the reduced order system. Experimental results are included using SRL’s 
second-generation, reduced-order spacecraft simulator with a state feedback linearized 
controller to demonstrate its ability to navigate a closed circular path with the proposed 
actuator configuration. A general overview of this spacecraft simulator is presented with 
additional details on the simulators being found in: Hall (2006), Eikenberry (2006), Price 
(2006), Romano and Hall (2006), Hall and Romano (2007a), and Hall and Romano 
(2007b). 

B. ROTO-TRANSLATION OF AN UNDERACTUATED SPACECRAFT 

With respect to the full order (6-DoF) system of roto-translation, it will be 
demonstrated that a set of two vectorable hemispherical thrusters in an opposing 
configuration can yield a minimum set of actuators for a controllable relative motion 
spacecraft system. The proposed actuator configuration can readily be seen to yield an 
underactuated mechanical system, that is to say the number of controls N u is fewer than 

the degrees of freedom. It will be shown that by adding a single-degree of freedom 
torque, the system can become fully actuated and a state-feedback linearizing controller 
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can be designed to drive the system to the desired position and orientation. Without the 
single-degree of freedom torque, this control actuator configuration equates to only two 
of the three body-fixed axes being directly torque controlled while the third can be 
indirectly controlled by capitalizing on the coupling of the terms of the moments of 
inertia for the system that appear in the Euler equations. Furthermore, given the cascaded 
nature of the angular motion equations whereby the angular rates drive the orientation, 
there has been much research to date to develop control algorithms to provide three-axis 
stabilization for an underactuated spacecraft with only two control torques. However, to 
date, there has not been a smooth time-invariant control algorithm to provide attitude 
stabilization of such a spacecraft with arbitrary inertia. By considering the attitude 
stabilization for the underactuated spacecraft system of angular motion equations in the 
general case, it may prove to be widely applicable to not only proximity operations but 
also de-tumbling and reorientation maneuvers of underactuated spacecraft that may either 
be designed with only two directly actuated control axes or be experiencing control 
actuator failure about one of its control axes during their mission life. 

C. STATE-OF-THE-ART IN UNDERACTUATED RIGID BODY 
STABILIZATION 

The problem of stabilization of a rigid spacecraft’s attitude dynamics and 
kinematics has been studied over the years in many papers and articles. However, the 
vast majority of the proposed control laws assume that the spacecraft is fully actuated. 
Wie and Barba (1985), Wie, Weiss and Arapostathis (1989), Vadali (1989), and Bajodah 
(2009a) address several nonlinear control techniques that provide time-invariant global 
asymptotic stability of the fully actuated spacecraft system of equations. Although these 
control laws provide for the necessary control of a nominally designed three-axis 
stabilized spacecraft in which three control torques exist, the question of control of 
underactuated spacecraft naturally enters when discussing actuator failures or when 
proposing minimally designed spacecraft systems. It is well understood that full order (3- 
DoF) control of the kinematics of such underactuated systems presents a challenging 
control problem; however, it should also be recognized that it has the distinct potential to 
provide several key benefits. Specifically, under the present thrust of Operationally 
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Responsive Space, where one focus is on transitioning from the typical large monolithic 
spacecraft design to one that embraces the new spacecraft paradigm of smaller, faster to 
manufacture and cheaper to produce and employ, underactuated control could provide a 
key enabling technology. Furthermore, in light of many current traditional spacecraft 
systems remaining in operation well past their intended design life despite actuator 
failures that degrade their capabilities, underactuated control could enable these aging 
systems to satisfy their original missions. 

The investigation of stabilization of underactuated spacecraft kinematics and 
attitude dynamics began with the theoretical establishment of the necessary and sufficient 
conditions for the controllability of a rigid body’s attitude with either gas thrusters or 
momentum exchange devices by Crouch (1984). He concluded that, for a spacecraft with 
momentum exchange devices, controllability is impossible with fewer than three, while 
for a spacecraft with independent paired jets, controllability is possible with two. It was 
later demonstrated by Kerai (1995), by using geometric control theory, that small-time 
local controllability of the rigid body equations assuming paired gas jets can indeed be 
achieved with only two control torques. Byrnes and Isidori (1991) proved that the full 
angular motion equations for a rigid spacecraft with only two controls cannot be 
asymptotically stabilized by smooth pure state feedback because they violate Brockett’s 
(1983) theorem on non-holonomic underactuated systems. With this in mind, they 
proposed a smooth feedback controller to affect partial stabilization of the rigid body 
model resulting in a revolute constant-rate motion about the uncontrolled axis of rotation. 
Later, Krishnan, McClamroch and Reyhanoglu (1995) proposed a hybrid control design 
combining continuous time features with discrete event features to affect a discontinuous 
feedback control strategy to stabilize any equilibrium attitude of an underactuated 
spacecraft with two momentum wheel actuators in finite time under the restriction that 
the total angular momentum vector of the system is zero. This control methodology 
translates directly to a study of an underactuated axi-symmetric spacecraft. Tsiotras, 
Corless, and Longuski (1995) and Tsiotras and Luo (2000) also dealt with control of 
underactuated axi-symmetric spacecraft by proposing a time-invariant feedback control 
law to asymptotically stabilize the orientation of two of the three body-fixed axes. In 
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addition to only providing for partial attitude stabilization of axi-symmetric spacecraft, 
their discussion was limited to cases where the angular velocity about the unactuated axis 
is zero at the start of the maneuver. Tsiotras and Schleicher (2000), and Tsiotras and 
Doumtchenko (2000) relaxed the restriction on the symmetry of the spacecraft slightly to 
consider a nearly axi-symmetric spacecraft by a small parameter and a set of time- 
invariant control laws are proposed to stabilize the angular velocity and attitude of a 
spacecraft about a certain axis by virtual control inputs of the two actuated angular rates. 

The global asymptotic rate stabilization problem without concern for kinematics 
of a fully asymmetric underactuated rigid spacecraft was addressed by Coverstone- 
Carroll (1996) through the use of a Variable Structure Controller (VSC). Bajodah 
(2009b) also addressed the rate-only stabilization problem for detumbling maneuvers 
through the use of singularly perturbed feedback linearization and generalized inverse 
control methodologies. Although both of these controllers prove to be robust to large 
initial angular velocities around all three axes in the presence of actuator torque 
limitations, they both require an additional controller to provide desired kinematic 
alignment after the detumbling maneuver. One such controller, as proposed by 
Coverstone-Carroll (1996), is a simple linear controller that is used to perform a series of 
eigenaxis rotations which precludes smooth attitude tracking. 

The problem of stabilization of both the kinematics and dynamics of an 
underactuated asymmetric spacecraft was most recently addressed by Casagrande, 
Astolfi, and Parisini (2008) who proposed a time-variant switching control law to effect 
global asymptotic stabilization of the closed-loop system. Although novel, the proposed 
law lacks detailed simulation results by considering only the case where the initial 
angular rates about two of the axes to include the unactuated axis are initially zero. 
Furthermore, real-world spacecraft with flexible parts, antennas, fuel slosh, etc., may 
preclude the use of time-variant control laws because they have the distinct potential of 
producing unacceptable transient response and might therefore lead to instability 
(Tsiotras and Doumtchenko, 2000). Behai et al. (2002) address the nonlinear tracking 
control of an axi-symmetric spacecraft by developing a kinematic controller to determine 
the desired actuated angular rates which are in turn used as control inputs to the dynamic 
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system through the use of standard back-stepping techniques. This method yields only 
asymptotic dynamic and kinematic stabilization results for an axi-symmetric rigid body 
given the restriction that the angular rate about the unactuated axis is close to zero but it 
does yield bounded results otherwise. 

The goal of this work is to extend the research into control algorithms for 
underactuated rigid spacecraft attitude control by proposing a novel, time-invariant 
smooth quaternion feedback regulator based on generalized inverse methods to affect 
three-axis attitude stabilization of the error quaternion kinematics for an underactuated 
rigid spacecraft with arbitrary inertia matrices and two realistically bounded body-fixed 
torques for required reorientation maneuvers. The problem of three-axis stabilization of 
the attitude of an underactuated spacecraft with direct control about only two of the body- 
fixed control axes will be addressed in the general case where the spacecraft’s attitude is 
referenced to an inertially fixed frame. From this, the proposed quaternion feedback 
regulator can be shown to be seamlessly modified to account for attitude control with 
respect to a relative motion frame of reference as is the case for a chaser-target 
rendezvous maneuver. After affecting attitude error stabilization, a spacecraft can be 
propelled towards another spacecraft via various navigation schemes such as 
conventional waypoint navigation. 

D. SCOPE OF THE DISSERTATION 

This dissertation advances the body of knowledge with respect to control of 
underactuated spacecraft in three key areas: 

1. Laboratory experimentation of a reduced-order underactuated spacecraft 
simulator with vectorable thrusters and a miniature control moment 
gyroscope. Using feedback linearizing control methodology coupled with 
Schmitt Trigger and Pulse Width Modulation logic, experimental results are 
presented which validate the capability of this novel control actuator design to 
propel the spacecraft simulator around a tightly constrained path. 

2. Analytical determination of the small-time local controllability of a generic 
full-order spacecraft under variations on the control inputs. This study is able 
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to provide the interested spacecraft systems engineer with the ability to 
determine the minimum number of control actuators necessary to maintain 
controllability. Furthermore, this analysis can be used to aid in dealing with 
both control actuator failures on existing spacecraft systems or planning for 
minimally designed spacecraft. 

3. A smooth time-invariant state feedback control logic based on quaternion 
feedback regulation is derived to yield stabilization of the error kinematics of 
a spacecraft with only two control torques and arbitrary inertia. Two separate 
control designs are presented, the first being Lyapunov function based and the 
other being perturbed feedback linearizing in nature. Results of the numerical 
results considering both of these designs are presented for various maneuvers. 
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II. LABORATORY EXPERIMENTATION OF GUIDANCE AND 
CONTROL OF SPACECRAFT DURING PROXIMITY 

MANEUVERS 

While presenting an overview of a robotic testbed for HIL experimentation of 
guidance and control algorithms for on-orbit proximity maneuvers, this chapter 
specifically focuses on exploring the feasibility, design and evaluation in a 3-DoF 
environment of a vectorable thruster configuration combined with optional miniature 
single gimbaled control moment gyro (MSGCMG) for an agile small spacecraft. 
Specifically, the main aims are to present and practically confirm the theoretical basis of 
small-time local controllability for this unique actuator configuration through both 
analytical and numerical simulations performed in previous works (Romano and Hall 
2006; Hall and Romano 2007a; Hall and Romano 2007b) and to validate the viability of 
using this minimal control actuator configuration on a small spacecraft in a practical way. 
Furthermore, the experimental work is used to confirm the controllability of this 
configuration along a fully constrained trajectory through the employment of a smooth 
feedback controller based on state feedback linearization and linear quadratic regulator 
techniques and proper state estimation methods. The chapter is structured as follows: 
First the design of the experimental testbed including the floating surface and the second 
generation 3-DoF spacecraft simulator is introduced. Then the dynamics model for the 
spacecraft simulator with vectorable thrusters and momentum exchange device are 
formulated. The controllability associated with this uniquely configured system is then 
addressed with a presentation of the minimum number of control inputs to ensure small 
time local controllability. Next, a formal development is presented for the state feedback 
linearized controller, state estimation methods, Schmitt trigger and Pulse Width 
Modulation scheme. Finally, experimental results are presented that demonstrate a 
closed-path circular trajectory about an arbitrary reference that is representative of a 
possible inspection of a target spacecraft by a given chaser. 
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A. 


THE NPS ROBOTIC SPACECRAFT SIMULATOR TESTBED 


Three generations of robotic spacecraft simulators have been developed at the 
NPS Spacecraft Robotics Laboratory, in order to provide for relatively low-cost HIL 
experimentation of GNC algorithms for spacecraft proximity maneuvers (see Figure 1). 
In particular, the second generation robotic spacecraft simulator testbed is used for the 
here-in presented research. The whole spacecraft simulator testbed consists of three 
components. The two components specifically dedicated to HIL experimentation in 3- 
DoF are a floating surface with an indoor pseudo-GPS (iGPS) measurement system and 
one 3-DoF autonomous spacecraft simulator. The third component of the spacecraft 
simulator testbed is a 6-DoF simulator stand-alone computer based spacecraft simulator 
and is separated from the HIL components. Additionally, an off-board desktop computer 
is used to support the 3-DoF spacecraft simulator by providing the capability to upload 
software, initiate experimental testing, receive logged data during testing and process the 
iGPS position coordinates. Figure 2 depicts the robotic spacecraft simulator in the 
Proximity Operations Simulator Facility (POSF) at NPS with key components identified. 
The main testbed systems are briefly described in the next sections with further details 
given in Hall (2006), Price (2006), Eikenberry (2006), Romano and Hall (2006), Hall and 
Romano (2007a), and Hall and Romano (2007b). 
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Figure 1. Three Generations of Spacecraft Simulators at the NPS Spacecraft 
Robotics Laboratory (first, second and third generations, from left to right) 


1. Floating Surface 

A 4.9 m by 4.3 m epoxy floor surface provides the base for the floatation of the 
spacecraft simulator. The use of planar air bearings on the simulator reduces the friction 
to a negligible level and with an average residual slope angle of approximately 2.6x10' 
deg for the floating surface, the average residual acceleration due to gravity is 
approximately 1.8x10' ms' . This value of acceleration is two orders of magnitude lower 
than the nominal amplitude of the measured acceleration differences found during 
reduced gravity phases of parabolic flights (Romano et al. 2007). 
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Figure 2. SRL’s 2nd Generation 3-DoF Spacecraft Simulator 


2. Three-DoF Robotic Spacecraft Simulator 

SRL’s second generation robotic spacecraft simulator is modularly constructed 
with three easily assembled sections dedicated to each primary subsystem. Prefabricated 
6105-T5 Aluminum fractional t-slotted extrusions form the cage of the vehicle while one 
square foot, 0.25 inch thick static dissipative rigid plastic sheets provide the upper and 
lower decks of each module. The use of these materials for the basic structural 
requirements provides a high strength to weight ratio and enable rapid assembly and 
reconfiguration. Table 1 reports the key parameters of the 3-DoF spacecraft simulator. 
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Subsystem 

Characteristic 

Parameter 

Structure 

Length and width 

0.30 m 


Height 

0.69 m 


Mass 

26 kg 


K 

0.40 kg-m 2 

Propulsion 

Propellant 

Compressed Air 


Equiv. storage capacity 

0.05 m 3 @ 31.03 MPa 


Operating pressure 

0.41 MPa 


Thrust (x2) 

0.159 N 


ISP 

34.3 s 


Total AF 

31.1 m/s 

Flotation 

Propellant 

Air 


Equiv. storage capacity 

0.05 m 3 @ 31.03 MPa 


Operating pressure 

.51 MPa 


Linear air bearing (x4) 

32 mm diameter 


Continuous operation 

-40 min 

CMG Attitude Control 

Max torque 

0.668 Nm 


Momentum storage 

0.098 Nms 

Electrical & Electronic 

Battery type 

Lithium-Ion 


Storage capacity 

12 Ah @ 28V 


Continuous Operation 

-6 h 


Computer 

1 PC 104 Pentium III 

Sensors 

Fiber optic gyro 

KVH Model DSP-3000 


Position sensor 

Metris iGPS 


Magnetometer 

MicroStrain 3DM-GX1 


Table 1. Key Parameters of the 2nd Generation 3-DoF Robotic Spacecraft 

Simulator 


a. Propulsion and Flotation Subsystems 

The lowest module houses the flotation and propulsion subsystems. The 
flotation subsystem is composed of four planar air bearings, an air filter assembly, dual 
4500 PSI (31.03 MPa) carbon-fiber spun air cylinders and a dual manifold pressure 
reducer to provide 75 PSI (0.51 MPa). This pressure with a volume flow rate for each air 
bearing of 3.33 slfm (3.33 x 10' m/min) is sufficient to keep the simulator in a friction- 
free state for nearly 40 minutes of continuous experimentation time. The propulsion 
subsystem is composed of dual vectorable supersonic on-off cold-gas thrusters and a 
separate dual carbon-fiber spun air cylinder and pressure reducer package regulated at 60 
PSI (0.41 MPa) and has the capability of providing the system 31.1 m/s A V . 
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b. Electronic and Power Distribution Subsystems 

The power distribution subsystem is composed of dual lithium-ion 
batteries wired in parallel to provide 28 volts for up to 12 Amp-Hours and is housed in 
the second deck of the simulator. A four port DC-DC converter distributes the requisite 
power for the system at 5, 12 or 24 volts DC. An attached cold plate provides heat 
transfer from the array to the power system mounting deck in the upper module. The 
current power requirements include a single PC-104 CPU stack, a wireless router, three 
motor controllers, three separate normally-closed solenoid valves for thruster and air 
bearing actuation, a fiber optic gyro, a magnetometer and a wireless server for 
transmission of the vehicle’s position via the pseudo-GPS system. 

c. Translation and Attitude Control System Actuators 

The 3-DoF robotic spacecraft simulator includes actuators to provide both 
translational control and attitude control. A full development of the controllability for 
this unique configuration of dual rotating thrusters and one-axis Miniature-Single 
Gimbaled Control Moment Gyro (MSGCMG) will be demonstrated in subsequent 
sections of this paper. The translational control is provided by two cold-gas on-off 
supersonic nozzle thrusters in a dual vectorable configuration. Each thruster is limited in 
a region ±7t/ 2 with respect to the face normal and, through experimental testing at the 
supplied pressure, has been demonstrated to have an ISP of 34.3 s and able to provide 
0.159 N of thrust with less than 10 msec actuation time (Lugini and Romano 2009). The 
MSGCMG is capable of providing 0.668 Nm of torque with a maximum angular 
momentum of 0.098 Nms. 

3. Six-DoF Computer-based Numerical Spacecraft Simulator 

A separate component of SRL’s spacecraft simulator testbed at NPS is a 6-DoF 
computer-based spacecraft simulator. This simulator enables full 6-DoF numerical 
simulations to be conducted with realistic orbital perturbations including aerodynamic, 
solar pressure and third-body effects, and earth oblateness up to J4. Similar to the 3-DoF 
robotic simulator, the numerical simulator is also modularly designed within a 
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MATLAB®/Simulink® architecture to allow near seamless integration and testing of 
developed guidance and control algorithms. Additionally, by using the 
MATLAB®/Simulink® architecture with the added Real Time Workshop™ toolbox, the 
developed control algorithms can be readily transitioned into C-code for direct 
deployment onto the 3-DoF robotic simulator’s onboard processor. A full discussion of 
the process by which this is accomplished and simplified for rapid real-time 
experimentation on the 3-DoF testbed for either the proprietary MATLAB® based 
XPCTarget™ operating system is given in (Hall 2006; Price 2006) or for an open-source 
Linux based operating system with the Real Time Application Interface (RTAI) is given 
in (Bevilacqua, Hall, Horning, and Romano 2009). 

B. DYNAMICS OF A 3-DOF SPACECRAFT SIMULATOR WITH 
VECTORABLE THRUSTERS AND MOMENTUM EXCHANGE DEVICE 

Two sets of coordinate frames are established for reference: an inertial frame I 

with orthogonal axes defined by the unit vectors j , i 2 , ( 3 j and body-fixed frame 2? with 

r /v /v /v 'j 

orthogonal axes defined by the unit vectors \b v b 2 ,b 3 1. These reference frames are 

depicted in Figure 4 along with the necessary external forces and parameters required to 
properly define the simulators motion. The origin of the body-fixed coordinate system is 
taken to be the center of mass C of the spacecraft simulator and is assumed to be 

collocated with the simulator’s geometric center. By the effects of the flat floor, b 3 is 

yv ^ 

maintained aligned with i 3 while b x is defined to be in line with the thrusters points of 
action. The position and velocity vectors of 2? with respect to I expressed in 2? are 
given by r and v so that r BI marks the position of the simulator with respect to the origin 

of I as measured by the inertial measurement sensors and provides the vehicle’s two 
degrees of translational freedom. The vehicle’s rotational freedom is described by an 

^ /V ^ 

angle of rotation 0 S between b x and i x about 6,. The angular velocity of 2? with respect 
to I expressed in 2? is thus limited to one degree of freedom and is denoted by co .. The 
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spacecraft simulator is assumed to be rigid and therefore a constant moment of inertia 
(7 3 ) exists about the 4. Furthermore, any changes to the mass of the simulator (m) due 
to thruster firing are neglected. 

The forces imparted at a distance L from the center of mass by the vectorable on- 
off thrusters are denoted by F a and F b respectively. The thrust angle a a defines the 

orientation of thruster a in 22 and is the angle measured from 4 in a clockwise 
direction (right-hand rotation) to F fl . Likewise, thrust angle a b defines the orientation of 
thruster b in $ and is the angle measured from -b x in a clockwise direction (right-hand 
rotation) to F /; . The torque imparted on the vehicle by a momentum exchange device 
such as the MSGCMG is denoted by T CMC and can be constrained to exist only about the 
4 axis. 



Figure 3. Schematic of SRL’s 2nd Generation Spacecraft Simulator 


1, Controllability Analysis of a Single-Gimbaled Control Moment 
Gyroscope for a Reduced Order System 

The feasibility and utility of using a single-gimbaled control moment gyroscope 
for the actuation of a 3-DoF spacecraft simulator can easily be demonstrated through 
development of the reduced order angular motion equations for the vehicle and the 
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MSGCMG. Beginning with the full order angular motion equations equipped with a 
momentum exchange device such as a MSGCMG given by Wie (1998, p. 437) 

^H + cox S H = T (1) 

where S H(?) e M 3 , V S H = [ s // p s H 2 , the total system’s angular momentum vector 
is expressed in the body-fixed control frame with respect to an inertial frame, 
Text- (0 e ® 3 ’ VT £AT =\_T E x T \'Text 2’Text i~\ ' s th e external torque vector acting on the 

body-fixed frame and co(t) e M 3 , Vco = [<y^<y 2 ,<y 3 ] r is the angular velocity vector of the 
body-fixed frame with respect to an inertial frame. 

The total system’s angular momentum vector includes both the rigid body’s 
angular momentum and the MSGCMG’s angular momentum. Therefore 

S H = J T o) + h (2) 

where J T e E 3a3 is the total moment of inertia of the spacecraft with the SGCMG such 
that 


J j J g l J 


CMC 


(3) 


where J B e M 3 * 3 


h 


J P + 


1 0 
0 1 
0 0 


0 

0 

2 

/77 r 

" l CMG'CMG 


(4) 


represents the sum of the second moment of inertia dyadic of the SGCMG platform 
J p e E 3j3 and the second moment of inertia of a point mass concentrated in the center of 
mass m CMG of the SGCMG with respect to the center of mass of the spacecraft r CMC 
assuming the SCCMG is aligned with the spacecraft’s third control axis. J CMC e E 3 * 3 is 
the second moment of inertia dyadic of the MSGCMG cluster represented in the body- 
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fixed frame. This represents the summation of the second moment of inertia dyadic of 
the MSGCMG body and includes both the contribution of its spinning inertia disc and 
associated hardware (Bevilacqua, Izzo, and Valente 2003). Finally, 

h(f) e M 3 , Vh = [/ip/ij./ij] 7 represents the total momentum vector of the SGCMG 
expressed in the body-fixed frame. Differentiation of (2) yields 


? H = J T d) + j T (i) + h 
= J T <a + of J T oj + h 


(5) 


where of e M 3 ' 3 represents the skew symmetric matrix with respect to co 


co 


x 


0 

-co 3 

co 2 

co 3 

0 

-co x 

-co 2 

co x 

0 


( 6 ) 


By combining (1), (2), and (5), we get 

/ 7 .cb + ft> x / 7 .co + h + cox(/ J .co + h) = T EXT (7) 

Furthermore, by introducing the control torque vector generated on the body-fixed 
frame by the MSGCMG denoted as T CMG (t) e M 3vl , VT CMC = \T CMGV T CMG2 ,T CMG2 ^ and a 

control vector T(0eU 3,1 ,VT = [T^T 2 ,T 3 f to represent the torque required by a given 
control law, (7) becomes 

J T d) + co' J T (o = T CMG + T EXT (8) 

where 

T CMG =-(<z> x / r (o + h + (oxh) = T (9) 

(9) can be solved for the time derivative of the angular momentum vector as 

h = -T-(o x / r (o -coxh (10) 
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v 

Figure 4. Orthonormal Vectors of MSGCMG Unit (After Kurokawa 1998) 

Figure 4 depicts the internal vectors of the MSGCMG unit. Three mutually 
orthogonal unit vectors exist where g is the gimbal vector, his the angular momentum 
vector, and c is the torque vector where 

3h , 

c = — = gxh (11) 

do 

The gimbal vectors are constant while the others are dependent upon the 
MSGCMG’s gimbal angle 8 . The total angular momentum becomes a general function 
of the CMG gimbal angles 5 and the constant angular momentum of the unit’s rotor 
wheel denoted by h w such that 



h = h w [0,cosc),sinc>] r (12) 

The total output torque of the MSGCMG in the absence of coupling terms due to 
spacecraft motion then can be found by taking the time derivative of (12) as 

T cmg = -h = [0, h w sin 88, -h w cos 88 J (13) 

By taking the third body-fixed axis to be the yaw axis as depicted in Figure 4, and 
considering only the components about this axis due to compensation of the other two 
axes by the reactions of the floor, (8) and (9) simplify to 

J 3 cb 3 -T 3 - -h w cos 88 (14) 
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To determine the steering logic for the instantaneous angular rate to command the 
MSGCMG’s gimbal motor, it is thus necessary to solve (14) for the gimbal rate 8 given 
the desired torque as generated by a given attitude control law. Once the gimbal angular 
position reaches ±x/2 , the MSGCMG cannot generate any torque about the third body- 
fixed axis and it is therefore in a saturated singular configuration. 

2. Reduced Order Dynamics with Dual Rotating Thrusters and 
MSGCMG 

The translation and attitude motion of the simulator are governed by the equations 

r «, = 

v = m~ x 'R b F 

\ Bi rn iV g A (15) 

0, = co 3 
C0 3 = J 3 % 

where B F(7) e Hr , V B F = B F n J represents the thrusters’ inputs limited to the region 

±7z/ 2 with respect to each face normal and T. (f) el is the attitude input. 1 R ls , 11 F and 
T 3 are given by 



c 0 3 

-S0 3 

s6*3 

c6 3 


(16) 


*F r = *FJ + *F s r = Kca,+F,ca„-F |S a. + F,s«,f (17) 

T 3 = [^cmg + ^ ( ~F a s ci a — F b s oc b )J (18) 

where s» = sin(»),c» = cos(*). 

The internal dynamics of the vectorable thrusters are assumed to be linear 
according to the following equations 


<=P u J a = KX,<=P h J h =J b X 


(19) 
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where J a and J b represent the moments of inertia about each thruster rotational axis 
respectively and T h eR represent the corresponding thruster rotation control 

input. 

The system’s state equation given by (15) can be rewritten in control-affine 
system form as (LaValle 2006) 

N„ 

x = f (x) + 2>,g, (x) = f (x) + G(x)u, x e W Nx (20) 

i=i 

where N u is the number of controls. With M :V ' representing a smooth N x -dimensional 
manifold defined be the size of the state-vector and the control vector to be in U :V " 
Defining the state vector xef 10 as \ T = [x p x 2 ,...,x 10 ] = 

10i/,i > Tv. 2 > > a u » a i, ’ v niA . v fl/,2 » °h > Pa > Ph I and the control vector ueU 5 as 

u r = [u v u 2 ,...,u 5 ] = [ B F a , B F b ,T CMG ,T a ,T b ], the system’s state equation becomes 


x f ( X ) + G(x)u [.* 6 ,x 7 ,x 8 ,x 9 ,x 10 ,0 1t5 ] + 


5.v5 


G,(x)_ 


u 


( 21 ) 


where the matrix G, (x) e R 5 ' 5 is obtained from (15) as 


G,(x) = 


-m l [c x 3 c ,v' 4 - s v 3 s x 4 ] 
-m 1 [ci 3 sx 4 +sj 3 cjc 4 ] 
- J 3 _1 s x 4 L 

0 

0 


m 

m 


1 [c x 3 c v 5 - s v 3 s a s ] 
1 [c a 3 s x 5 + s x 3 c x 5 ] 
- 7 3 1 s x 5 L 

0 

0 


0 

0 

0 


0 

0 

0 

J-' 


0 

0 

0 

0 


0 0 j; 


( 22 ) 


With the system in the form of (20) given the vector fields in (21) and (22), and 
given that f(x) e M 1() (the drift term) and G(x)el 10 ' 5 (the control matrix of control 
vector fields) are smooth functions, it is important to note that it is not necessarily 
possible to obtain zero velocity due to the influence of the drift term. This fact places the 
system in the unique subset of control-affine systems with drift and, as seen later, will 
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call for an additional requirement for determining the controllability of the system. 
Furthermore, when studying controllability of systems, the literature to date restricts the 
consideration to cases where the control is proper. Having a proper control implies that 
the affine hull of the control space is equal to U :V “ or that the smallest subspace of U is 
equal to the number of control vectors and that it is closed (Sussman 1987; Sussman 
1990; Bullo and Lewis 2005; LaValle 2006). With a system such as a spacecraft in 
general or the simplified model of the 3-DoF simulator in particular, the use of on-off 
cold-gas thrusters restrict the control space to only positive space with respect to both 
thrust vectors leading to an unclosed set and thus improper control space. In order to 
overcome this issue, a method which leverages the symmetry of the system is used by 
which the controllability of the system is studied by considering only one virtual rotating 
thruster that is positioned a distance L from the center of mass with the vectored thrust 
resolved into a y and x-component. In considering this system perspective, the thruster 
combination now spans M 2 and therefore is proper and is analogous to the planar body 
with variable-direction force vector considered in (Lewis and Murray 1997; Bullo and 
Lewis 2005). Furthermore, under the assumption that the control bandwidth of the 
thrusters’ rotation is much larger than the control bandwidth of the system dynamics, the 
internal dynamics of the vectorable thrusters can be decoupled from the state and control 
vectors for the system yielding a thrust vector dependent on simply a commanded angle. 
Thus the system’s state vector, assuming that both thrusters and a momentum exchange 
device are available, becomes x r = [jq,jc 2 ,...,x 6 ] = [r BI v r BI 2 ,# 3 ,v B/ ,,v B/ 2 ,&> 3 ] e R 6 and 

the control vector is u T =[u l ,u 2 ,u 3 ] = [ B F 1 , B F 2 ,T 3 \ = U i so that the system’s state 
equation becomes 


x = f (x) + G(x)u = [.r 4 ,x 5 ,x 6 ,0,0,0] r + 


'3.v3 


G,(x)_ 


u 


(23) 


where f(x) e M 6 and G(x) e M 6a3 are again smooth functions. The matrix G 1 (x)can be 
obtained by considering the relation of the desired control vector to the body centered 
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reference system, in the two cases of positive force needed in b x ( B U l > 0 ) and negative 
force needed in b , ( R U ] <0). In this manner, the variables in (22) and (23) can be 
defined as 


B t/, < 0 - 
B U l > 0 ■ 


[u 7 = [ B F l , b F 2 ,T 3 ] = [~F a c x 4 ,-F a s x 4 , T cmg ] 
[d = -L,F b = 0 

[u 7 =[ B F l , B F 2 ,T 3 \ = [F b cx 5 ,F b sx 5 ,T CMG ] 

\d = L,F= 0 


(24) 


yielding the matrix in 


(x)' 


D 3*3 


through substitution into (22) as 


G.(x) = 


m c x 3 -m s x 3 


0 


m 1 s x 3 

0 


m cx 3 0 
-dJ, 1 /(' 


(25) 


When the desired control input to the system along b x is zero, both thrusters can 
be used to provide a control force along b 2 , while a momentum exchange device provides 
any required torque. In this case, the control vector u becomes 
u r = [u v u 2 \ = [ b F 2 ,T 3 ] = U 2 and control matrix of control fields becomes G(x) e M 6x2 in 
(23) such that the variables in (22) and (23) can be defined as 


h U 2 =0^{ 


u r =[ B F 2 ,T 3 ] = [Fsa,T CMG \ 

F = F a =F b ,cc = x A =-x 5 =--sign( B U 2 ) 


(26) 


yielding the matrix G, (x) e R 3x ° through substitution into (22) as 


Gj(x) 


-2m 1 s x 3 0 
2 m~ l cx 3 0 

0 7 3 1 


(27) 
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As will be demonstrated later, the momentum exchange device is not necessary to 
ensure small time controllability for this system. In considering this situation, which also 
occurs when a control moment gyroscope is present but is near the singular conditions 
and therefore requires desaturation, the thruster not being used for translation control can 
be slewed to ±nj2 depending on the required torque compensation and fired to affect the 
desired angular rate change. The desired control input to the system with respect to 
b x ( B U l ) can again be used to define the desired variables such that 


B U X < 0 -» 
B U, > 0 -» 


ju T =[ B F l , B F 2 ,T 3 ] = [ B F 1 , b F 2 , B F 2lrot3 ] = [~F a c x 4 ,-F a sx 4 , F h d sx 5 ] 
[d = -L, x 5 = + ;r/2 

ju r = [ ■ b F x , b F 2 , T 3 ] = [ B F l , b F 2 , B F 2lrot3 ] = [F b c x 5 , F b s x 5 , -F a d s x 4 ] 

[d = L, x 4 = ±7r/ 2 


(28) 


yielding the matrix G, (x) e R 3r3 through substitution into (22) as 


G,(x) 


m 

cx 3 

—m 

sx 3 

nf 

1 sx 3 

m 1 

cx 3 


0 -dJ , 1 


- ( md) sx 3 
(nid) 'cx 3 



(29) 


In case of zero force requested along x with only thrusters acting, the system 
cannot in general provide the requested torque value. 

A key design consideration with this type of control actuator configuration is that 
with only the use of an on/off rotating thruster to provide the necessary torque 
compensation, fine pointing can be difficult and more fuel is required to affect a desired 
maneuver involving both translation and rotation. 
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C. SMALL-TIME LOCAL CONTROLLABILITY 

Before studying the controllability for a nonlinear control-affine system of the 
form in (20), it is important to review several definitions. First, the set of states reachable 
in time at most T is given by R z ( T) by solutions of the nonlinear control-affine 
system. 

Definition 1 (Accessibility) 

A system is accessible from x 0 (the initial state) if there exists r > 0 such that the 
interior of R z (x 0 ,< t ) is not an empty set for t e ]0, r] (Bullo and Lewis 2005). 

Definition 2 (Proper Small-time Local Controllability) 

A system is small-time locally controllable (STLC) from x 0 if there exists r > 0 
such that x 0 lies in the interior of /? x (x 0 ,< t) for each t e ]0, r]for every proper control 
set U (Bullo & Lewis, 2005). Assuming that at x(0) = 0 this can also be seen under time 
reversal as the equilibrium for the system x () can be reached from a neighborhood in 
small time (Sussman 1987; Sussman 1990). 

Definition 3 (Proper Control Set) 

A control set u 7 =[u v ...,u k ] is termed to be proper if the set satisfies a constraint 

uef where affinely spans LJ 7 (Sussman 1990; Bullo and Lewis 2005; LaValle 
2006). 

Definition 4 (Lie derivative) 

The Lie derivative of a smooth scalar function g(x) e R with respect to a smooth 
vector field f (x) e R' v ' is a scalar function defined as (Slotine 1991, p. 229) 


L { g = Vg f 


r 


/i(x) 

dg 

dg 

dx. 

dx N 


1 

ly x 

A (x) _ 


(30) 
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Definition 5 (Lie Bracket) 


The Lie bracket of two vector fields f (x) e and g(x) e Wl N] is a third vector 
field [f,g] e defined by [f, g] = Vg f-Vf g, where the i-th component can be 
expressed as (Slotine 1991, p. 230) 


[ f ’g],= Z 


N ' f f . 

J8x J 1 dx u 


i =i 


(31) 


Using Lie bracketing methods that produce motions in directions that do not seem 
to be allowed by the system distribution, sufficient conditions can be met to determine a 
system’s STLC even in the presence of a drift vector as in the equations of motion 
developed above. These sufficient conditions involve the Lie Algebra Rank Condition 
(LARC). 

Definition 6 (Associated Distribution A(x)) 

Given a system as in (6), the associated distribution A(x) is defined as the vector 
space (subspace of ) spanned by the system vector fields Cg^.-g^ . 

Definition 7 

The Lie algebra of the associated distribution .Z’ (A) is defined to be the 

distribution of all independent vector fields that can be obtained by applying subsequent 
Lie bracket operations to the system vector fields. Of note, no more than N x vector 

fields can be produced (LaValle 2006). With dim[X(A)] < N x , the computation of the 

elements of ./’(A) ends either when N T independent vector fields are obtained or when 
all subsequent Lie brackets are vector fields of zeros. 

Definition 8 (Lie Algebra Rank Condition [LARC]) 

The Lie Algebra Rank Condition is satisfied at a state x if the rank of the matrix 
obtained by concatenating the vector fields of the Lie algebra distribution at x is equal to 
N x (the number of state). 
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For a driftless control-affine system, following the Chow-Rashevskii Theorem, 
the system is STLC if the LARC is satisfied (Lewis and Murray 1997; Bullo and Lewis 
2005; LaValle 2006). However, given a system with drift, in order to determine the 
STLC, the satisfaction of the LARC it is not sufficient: In addition to the LARC, it is 
necessary to examine the combinations of the vectors used to compose the Lie brackets of 
the Lie algebra. From Sussman’s General Theorem on Controllability, if the LARC is 
satisfied and if there are no ill formed brackets inX(A), then the system is STLC from 

its equilibrium point (Sussman 1987). The Sussman’s theorem, formally stated is 
reported here below. 

Theorem 1 (Sussman’s General Theorem on Controllability) 

Consider a system given by (20) and an equilibrium point p e M /V,Jl such that 
f (p) = 0 . Assume £ (A) satisfies the LARC at p . Furthermore, assume that whenever 
a potential Lie bracket consists of the drift vector f (x) appearing an odd number of 
times while g,(x),...^^ (x) all appear an even number of times to include zero times 

(indicating an ill formed Lie bracket), there are sufficient successive Lie brackets to 
overcome this ill formed Lie bracket to maintain LARC. Then the system is STLC from 
p (Sussman 1987; Sussman 1990). 


As it is common in literature, an ill formed bracket is dubbed a “bad” bracket 
(Sussman 1987; Sussman 1990; Lewis and Murray, 1997, Bullo and Lewis, 2005; 
LaValle 2006). Conversely, if a bracket is not “bad,” it is termed “good”. As an 
example, for a system with a drift vector and two control vectors, the bracket 
[f, [gj, gj ]] is bad, as the drift vector occurs only once while the first control vector 


appears twice and the second control vector appears zero times. Similarly, the bracket 




is good as the first control vector appears only once. Therefore, it can be 


summarized that if the rank of the Lie algebra of a control-affine system with drift is 
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equal to the number of states and there exist sufficient “good” brackets to overcome the 
“bad” brackets to reach the required LARC rank, then the system is small time locally 
controllable. 

D. SMALL-TIME LOCAL CONTROLLABILITY CONSIDERATIONS FOR 

THE 3-DOF SPACECRAFT SIMULATOR 

The concept of small time local controllability is better suitable than the one of 
accessibility for the problem of spacecraft rendezvous and docking, as a spacecraft is 
required to move in any directions in a small interval of time dependent on the control 
actuator capabilities (e.g., to avoid obstacles). The finite time r can be arbitrary if the 
control input is taken to be unbounded and proper (Sussman 1990; Bullo and Lewis 2005; 
LaValle 2006). 

While no theory yet exists for the study of the general controllability for a non¬ 
linear system, the STLC from an equilibrium condition can be studied by employing 
Sussman’s theorem. For the case of spacecraft motion, in order to apply Sussman’s 
theorem, we hypothesize that the spacecraft is moving from an initial condition with 
velocity close to zero (relative to the origin of an orbiting reference frame). 

In applying Sussman’s General Theorem on Controllability to the reduced system 
equations of motion presented in (23) with G, (x) given in (25), the Lie algebra evaluates 
to 


£ (A) = span {g,, g 2 , g 3 [f, g, ] , [f, g 2 ], [f, g 3 ]} (32) 

so that dim[^(A)] = N x = 6. In order to verify that this is the minimum number of 

actuators required to ensure STLC, the Lie algebra is reinvestigated for each possible 
combination of controls. Appendix A includes a developed MATLAB® function for 
determining the Lie algebra for an arbitrary nonlinear system. The resulting analysis, as 
summarized in Table 2, demonstrates that the system is STLC from the systems 
equilibrium point at f (p) = 0 given either two rotating thrusters in complementary semi¬ 
circle planes or fixed thrusters on opposing faces providing a normal force vector to the 
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face in opposing directions and a momentum exchange device about the center of mass. 
For instance, in considering the case of control inputs u 7 =[m p m 2 ] = [ b F v ,T,] e U 2 , (23) 
given (27) becomes 

X = f(x) + g,(x)M 1 +g 2 (x)M 2 

r nr r -l -l ~\ T r -i“l r ^3) 

= [v 4 ,x 5 ,.x; 6 ,0,0,0j +[^0,0,0,2m sx 3 ,2m cx 3 ,0j u t + ^0,0,0,0,0, J 3 J u 2 


with the equilibrium point p such that f(p) = 0 is p = [x p jc 2 ,jc 3 ,0,0,0] T . The ./’(A) 
formed by considering the associated distribution A(x) and successive Lie brackets as 


f, 

go 

[ f ’gl]’ 

[gpg 2 ]> 

[ f ’[ f ’gl]]’ 

[f,[gl,g 2 ]]. 

[gp[f>gl]]> 

[gl.[gpg 2 ]]. 

[g 2 ,[f,gl]]. 

[g 2 ,[gl,g 2 ]]. 

[f,[f,[f,g 2 ]]], 

[f,[f,[g„g 2 ]]], 


g 2 

[ f ’g 2 ] 

[f-[f ,g 2 ]], 

[gi-[ f .g 2 ]]> 

[g 2 ,[f,g 2 ]], [f,[f,[f,gi]]], 

[f,[f,[f,g 2 ]]], [f,[g P [f,g 2 ]]] 


The sequence can first be reduced by considering any “bad” brackets in which the 
drift vector appears an odd number of times and the control vector fields each appear an 
even number of times to include zero. In this manner the Lie brackets [gi,[f,gj]] 

and [g 2 , [f, g 2 ]] can be disregarded. 

By evaluating each remaining Lie bracket at the equilibrium point p, the linearly 
independent vector fields can be found as 
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§1 = [® 3 xl’~ m lsx 3 > m 1 CX 3 ,-J 3 1 L^ 

§2 = [®5jc1’^3 ] 

[ f . g t ] = Vg, ■ f - Vf ’ gi = \_m~ l sx 3 ,-m~ 1 cx 3 , 0 4xl J 

[f, g 2 ] = V §2 • f - Vf ' §2 = [02,1, -J 3 1 » °3,, J 

[gp [ f . g 2 ]] = v [ f . g 2 ] • gi - Vgi ■ [f,g 2 ] = [o M , ~{rnJ 3 )“' cx 3 , - ( mJ 3 )"' sx 3 ,0] 

[f, [gp [f, g 2 ]]] = v [gi, [ f ,g 2 ]]• f - Vf • [g,, [f ,g 2 ]] = [2 L(mJ 3 ) _1 cx 3 ,2L( mJ 3 )~‘ s* 3 ,0 4vl 

(34) 


Therefore, the Lie algebra comprised of these vector fields is 


^’(A) = 5pan{g 1 ,g 2 ,[f,g 1 ],[f,g 2 ],[g l ,[f,g 2 ]],[f,[g 1 ,[f,g 2 ]] 


(35) 


yielding dim [./"(A)] = N x = 6, and therefore the system is small time locally 
controllable. 


Control Inputs 

Thruster Positions 

dim(.£(A)) 

Controllability 

B F i 

a a =a h = 0 

2 

Inaccessible 

B F 2 

««=-«* =+V 2 

5 

Inaccessible 

T 

1 CMG 

NA 

2 

Inaccessible 

b f t 

1 1 ’ ± CMG 

a a = a h = 0 

6 

STLC 

( b F T \or( B F b F ) 

\ 1 2’ ± CMGJ ur \ 1 2’ 1 2/rot3) 

=+^/ 2 

6 

STLC 

B f B J7 

1 1» 1 2 

\a a \<7r/2,\a h \<7r/2 

6 

STLC 


Table 2. STLC Analysis for the 3-DoF Spacecraft Simulator 

E. NAVIGATION AND CONTROL OF THE 3-DOF SPACECRAFT 
SIMULATOR 


In the current research, the assumption is made that the spacecraft simulator is 
maneuvering in the proximity of an attitude stabilized target spacecraft and that this 
spacecraft follows a Keplarian orbit. Furthermore, the proximity navigation maneuvers 
are considered to be fast with respect to the orbital period. A pseudo-GPS inertial 
measurement system by Metris, Inc. (iGPS) is used to fix the ICS in the laboratory setting 
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for the development of the state estimation algorithm and control commands. The X-axis 
is taken to be the vector between the two iGPS transmitters with the Y and Z axes forming 
a right triad through the origin of a reference system located at the closest corner of the 
epoxy floor to the first iGPS transmitter. Navigation is provided by fusing of the 
magnetometer data and fiber optic gyro through a discrete Kalman filter to provide 
attitude estimation and through the use of a linear quadratic estimator to estimate the 
translation velocities given inertial position measurements. Control is accomplished 
through the combination of a state feedback linearized controller, a linear quadratic 
regulator, Schmitt trigger logic and Pulse Width Modulation using the minimal control 
actuator configuration of the 3-DoF spacecraft simulator. Figure 5 reports a block 
diagram representation of the control system. 



Figure 5. Block diagram of the control system of the 3-DoF spacecraft simulator, 

(vi, v 2 , v 3 control inputs to the system) 
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1. Navigation Using Inertial Measurements with Kalman Filter and 
Linear Quadratic Estimator 

In the presence of the high accuracy, low noise, high bandwidth iGPS sensor with 
position accuracy to within 5.4 mm with a standard deviation of 3.6 mm and 
asynchronous measurement availability with a nominal frequency of 40 Hz, a full-order 
linear quadratic estimator with respect to the translation states is implemented to 
demonstrate the capability to estimate the inertial velocities in the absence of 
accelerometers. Additionally, due to the effect of noise and drift rate in the fiber-optic 
gyro, a discrete-time linear Kalman filter is employed to fuse the data from the 
magnetometer and the gyro. Both the gyro and magnetometer are capable of providing 
new measurements asynchronously at 100 Hz. Experimental results are presented in 
subsequent sections to demonstrate the filter’s effectiveness. 


a. Attitude Discrete-time Kalman Filter 

With the attitude rate being directly measured, the measurement process 
can be modeled in state-space equation form as: 







(36) 



where a) is the measured gyro rate, /? is the gyro drift rate, rj and Tj Pg are the 
associated gyro output measurement noise and the drift rate noise respectively. y/ m is the 
measured angle from the magnetometer, and is the associated magnetometer output 
measurement noise. It is assumed that n, and /? are zero-mean Gaussian white- 
noise processes with variances given by cr^, al i; and cr m respectively. Introducing 
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the state variables x 7 control variables u = co g , and error variables 

w r = ,}/ fjg J and v = ?j l/m , (36) and (37) can be expressed compactly in matrix form as 

x(t) = A(t)x(t) + fi(t)u(t) + G(t)w(t) (38) 

z (t) =//x(t) + v(t) (39) 

In assuming a constant sampling interval At in the gyro output, the system 
equation (38) and observation equations (39) can be discretized and rewritten as 

X * + l =®A +r >A :+ T i.W, (40) 

z * = tf* x * +v * ( 41 ) 


where 




= £' AA/ 


1 -At 
0 1 


(42) 


and 


r 





(43) 


The process noise covariance matrix used in the propagation of the 
estimation error covariance given by (Gelb 1974; Crassidis and Junkins 2004) 

f = J^ 1 J f 4+1 0(t k+l ,T)G(r)E [w{r)w T {a)]G T {a)(b{t k+v a)drda (44) 


can be properly numerically estimated given a sufficiently small sampling interval by 
following the numerical solution by van Loan (Crassidis and Junkins 2004). First, the 
following 2/7 x 2/7 matrix is formed: 


-A GQG t 

0 a t 


At 


(45) 
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where At is the constant sampling interval, A and G are the constant continuous-time 
state matrix and error distribution matrix given in (38), and Q is the constant continuous¬ 
time process noise covariance matrix 


£> = £'{w(0w r (0} 




The matrix exponential of (45) is then computed by 


X 

^ 2 "| 


X 


0 

%_ 


0 

_ 


(46) 


(47) 


where is the state transition matrix from (42) and Q k Y k )' ■ Therefore, the 

discrete-time process noise covariance is 


^=( m ^) r =^ 2 


l/3a^At 3 +a> -l/2a^Ar 
-\Jla 2 fig At 2 <j- p gA t 


(48) 


The discrete-time measurement noise covariance is 

r k = E {v k vl} = (J L ( 49 ) 

Given the filter model as expressed in (36) and (37), the estimated states 
and error covariance are initialized where this initial error covariance is given 
byP 0 = £’|x(t 0 )x r (t 0 )|. If a measurement is given at the initial time, then the state and 
covariance are updated using the Kalman gain formula 

+ (50) 


where P' k is the a priori error covariance matrix and is equal to P 0 . The updated or a 
posteriori estimates are determined by 
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( 51 ) 


r k = x k + K k [zj. - H k x k ] 
K=[h x i~K k H k ]P; 


where again with a measurement given at the initial time, the a priori state x k is equal to 
x 0 . The state estimate and covariance are propagated to the next time step using 


x * + i=®A+ r A 


(52) 


If a measurement is not given at the initial time step, or any time step 
during the process, the estimate and covariance are propagated to the next available 
measurement point using (52). 


b. Translation Linear Quadratic Estimator 

With the measured translation state from the iGPS sensor being given by 


z = 


10 0 0 

0 1 0 0 

c 



(53) 


the dynamics of a full-order state estimator is described by the equation 

x = x - x = [ Ax + 5u] - [~Ax + Bu + L lqe (z - Cx)] 

= A(x-x)-L lqe (Cx-Cx) (54) 

— (^A — L L q E C^x 


where 


Ax + Bu : linearized plant dynamics 

Ax + Bu : system model 

L lqe : linear quadratic estimator gain matrix 

Cx : predicted measurement 
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The observer gain matrix L lqe c an be solved using standard linear 
quadratic estimator methods as (Bryson 1993) 

L lqe =PC t R t 1 (55) 

where P is the solution to the algebraic Riccati equation 

AP + PA t - PC T R r l CP + Q t = 0 (56) 

and 2 r and R r arc the associated weighting matrices with respect to the translational 
degree of freedom defined as 

Q t = diag (1 / ||Ar max f , 1 / 1 | Ar max f , 1 /1|Av max f , 1 / 1 | Av max f) 

(57) 

R I =diag(ll(F_J 2 M(F lm f) 

where ||Arl|,||Avl| are taken to be the maximum allowed errors between the current 
and estimated translational states and F max is the maximum possible imparted force from 
the thrusters. 

Table 3 lists the values of the attitude Kalman filter and translation state 
observer used for the experimental tests. 
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2. Smooth Feedback Control via State Feedback Linearization and 
Linear Quadratic Regulation 

Considering a Multi-Input Multi-Output (MIMO) nonlinear system in control- 
affine form, the state feedback linearization problem of nonlinear systems can be stated 
as follows: obtain a proper state transformation 

z = ®(x) where z e M. N * (58) 

and a static feedback control law 

u = a(x)+/?(x)v where vel A, “ (59) 

such that the closed-loop system in the new coordinates and controls become 



is both linear and controllable. The necessary conditions for a MIMO system to be 
considered for input-output linearization are that the system must be square or N u = N y 
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where N u is defined as above to be the number of control inputs and N y is the number 
of outputs for a system of the expanded form (Isidori 1989; Slotine 1990) 


x = f (x) + G(x)u 


y = XM x ) =h(x > 

i=\ 


(61) 


The input-output linearization is determined by differentiating the outputs y t in 
(61) until the inputs appear. Following the method outlined in Slotine (1990, chap. 6) by 
which the assumption is made that the partial relative degree r is the smallest integer 

such that at least one of the inputs appears in y >'', then 


y t M = L r r % + X L zj L f r, - lh i ( x ) u j ( 62 ) 

7=1 

with the restriction that L g (x) ^ 0 for at least one j in a neighborhood of the 

equilibrium point x () . Letting 


£(*) = 


4,4' H 4W 
4i 4”"'4(x) 


4,4' H 4W 
4 , 4 ' 4 ( x ) 


4,4'"’X M 


L L, 
f 


r N _. -1 


k N y ( X ) 


(63) 


so that (63) is in the form 


1 

_1 


V*.W 

(r 2 ) 

y 2 ; 

— 

VM X ) 

(%,) 

[_ y N y ’] 


i 

_i 


+ 


£(x)u 


(64) 


the decoupling control law can be found where the N y xN y matrix E(\) is invertible 

over the finite neighborhood of the equilibrium point for the system as 
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( 65 ) 


v, -l/AO) 
v 2 - L t ' 2 h 2 (\) 

v N y ~L f rNy h Ny (x) 

With the above stated equations for the simulator dynamics in (23) given G, (x) 
as defined in (25), if we choose 

( 66 ) 

the state transformation can be chosen as 

z T = [\ (x), (x), /ij (x), Z^/ij (x), (x), L f h 3 (x)] = [r B/ ,0 3 ,\ BI ,co 3 ] (67) 

where z(t) e M 6 , Vz = [z,,z 2 ,...,z 6 ] r are new state variables, and the system in (23) is 
transformed into 

z = [z 4 , z 5 , z 6 , m 1 (c z 3 B F l - s z 3 B F 2 ), m l (s z 3 B F 1 + c z 3 B F 2 ), J 3 % ] (68) 

The dynamics given by (23) considering the switching logic described in (24), 
(26) and (28) can now be transformed using (68) and the state feedback control law 

[*F,r : ] = E(xf(v-b) (69) 

into a linear system 




(70) 


where 


b= L'!/?, (x),Lf 2 /z ? (x),Lf 3 /z 3 (x)^j 


(71) 
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and £(x) given by (63) with equivalent inputs v(t) e M 3 , Vv = [v 1 ,v 2 ,v 3 ] / and relative 

degree of the system at the equilibrium point x () is (r p r 2 ,r 3 ) = (2,2,2) . Therefore the 

total relative degree of the system at the equilibrium point, which is defined as the sum of 
the relative degree of the system, is six. Given that the total relative degree of the system 
is equal to the number of states, the nonlinear system can be exactly linearized by state 
feedback and with the equivalent inputs v ( , both stabilization and tracking can be 

achieved for the system without concern for the stability of the internal dynamics (Slotine 
1990). 

One of the noted limitations of a feedback linearized based control system is the 
reliance on a fully measured state vector (Slotine 1990). This limitation can be overcome 
through the employment of proper state estimation. HIL experimentation on SRL’s 
second generation robotic spacecraft simulator using these navigation algorithms 
combined with the state feedback linearized controller as described above coupled with a 
linear quadratic regulator to ensure the poles of (70) lie in the open left half plane 
demonstrate satisfactory results as reported in the following section. 

a. Feedback Linearized Control Law with CMG Rotational Control 
and Thruster Translational Control 

By applying (69) to the dynamics in (23) given G, (x) as defined in (25) 
where the system is taken to be observable in the state vector y = \r BI ,0 3 f = [x 3 ,x 2 ,x 3 ] r 
and by using thruster b for translational control (i.e., for the case B U l > 0 where 
B U l = v, c d 3 + v 2 s 0 3 and B U 2 = -v, s0 3 + v 2 c0 3 ), the feedback linearized control law is 

u T = [ B F l , b F 2 , T 3 ] = [m B U X , m B U 2 , rnL B U 2 + J 3 v 3 ] (72) 

which is valid for all x in a neighborhood of the equilibrium point x 0 . Similarly, the 
feedback linearized control law when B U 1 < 0 (thruster a is providing translation control) 

u r = [ B F l , b F 2 , T 3 ] = [m B U l , m B U 2 , -mL B U , + J 3 v 3 ] (73) 
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Finally, when B U 3 = 0 (both thrusters used for translational control) given 
G, (x) as defined in (27) is 

u t =[ b F 2 ,T 3 ] = [m b U 2 12 , J 3 v 3 ] (74) 

b. Feedback Linearized Control Law for Thruster Roto- 
Translational Control 

As mentioned previously, by considering a momentum exchange device 
for rotational control, momentum storage must be managed. For a control moment 
gyroscope based moment exchange device, desaturation is necessary near gimbal angles 
of 7i/ 2. In this region, due to the mathematical singularity that exists, very little torque 
can be exchanged with the vehicle and thus it is essentially ineffective as an actuator. To 
accommodate these regions of desaturation, logic can be easily employed to define 
controller modes as follows: If the MSGCMG is being used as a control input and if the 
gimbal angle of the MSGCMG is greater than 75 degrees, the controller mode is switched 
from normal operation mode to desaturation mode and the gimbal angle rate is directly 
commanded to bring the gimbal angle to a zero degree nominal position while the 
thruster not being directly used for translational control is slewed as appropriate to 
provide torque compensation. In these situations, the feedback linearizing control law for 
the system dynamics in (23) given G ] (x) as defined in (29) where thruster b is providing 

translational control ( 11 U . > 0), and thruster a is providing the requisite torque is 

u r = [ B F 1 , b F 2 , T 3 ] = \m B U l ,{mL B U 2 - J 3 v 3 ) /lL ,(mL B U 2 + J 3 v 3 ) /l\ (75) 

Similarly, the feedback linearizing control law for the system assuming 
thruster a is providing translational control ( B U l < 0) while thruster b provides the 
requisite torque is 

u r = [ B F 1 , b F 2 , T 3 ] = \m B U { , (mL B U 2 + / 3 v 3 )/lL , (mL B U 2 - J 3 v 3 ) /l\ (76) 
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3. Determination of the Thruster Angles, Forces and CMG Gimbal 
Rates 

In either mode of operation, the pertinent decoupling control laws are used to 
determine the commanded angle for the thrusters and whether or not to open or close the 
solenoid for the thruster. For example, if B U X > 0 , (72) or (75) can be used to determine 
the angle to command thruster b as 

a h ~ tan -1 ( b F 2 / B F l ) (77) 


and the requisite thrust as 

F b = yl B F 1 2 + B F 2 2 (78) 

If the MSGCMG is being used, the requisite torque commanded to the CMG is 
taken directly from (72). In the normal operation mode, with the commanded angle for 
thruster one not pertinent, it can be commanded to zero without affecting control of the 
system. Similarly, if B U i <0, (73) or (76) can be used to determine the angle to 

command thruster b and the requisite thrust analogous to (77) and (78). The requisite 
torque commanded to the CMG is similarly taken directly from (73). The required CMG 
torques can be used to determine the gimbal rate S CMG to command the MSGCMG by 
solving (14) given a desired body-fixed torque T 3 such that 

4 MG = ~ T i/(K cos 8 cmg ) (79) 

where h w is the constant angular momentum of the rotor wheel and 8 CMG is the current 
angular displacement of the wheel’s rotational axis with respect to the horizontal. 

If the momentum exchange device is no longer available and B U X > 0 , the thruster 
angle commands and required thrust value for the opposing thruster can be determined by 
using (75) as 

a a =-n/2sign(mL B U 2 + J 3 v 3 \ (80) 
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and 


F ~ -sign{a a )[mL B U 2 + J 3 v 3 )/L (81) 

given F 3 = —F a sa a L . Likewise, the thruster angle commands and required thrust value 
for the opposing thruster given B U 2 < 0 can be determined by using (76) as 

a b = 7r/2sign(mL B U 2 - J 3 v 3 \ (82) 

and 

F b = sign{cc b ){mL B U 2 -7 3 v 3 )/ L (83) 

given F 3 = -F b sin a b L . 

4. Linear Quadratic Regulator Design 

In order to determine the linear feedback gains used to compute the requisite 
equivalent inputs v ; to regulate the three degrees of freedom, so that 

lim z, (t) = r BI l (t) = r lim z 2 (t) = r B: (t) = r ref , 

t —>00 7 J t —^oo 7 J 

lim z 4 (t) = v BIA (; t ) = v ref , lim z 5 (t) = v B 2 (t) = v ref , 

t —>00 7 J t ->00 7 J 

a standard linear quadratic regulator is employed where the state-feedback law 
v = -Kz minimizes the quadratic cost function 

oo 

J (\) = ^{z T Qz + \ T Ry^dt (85) 

0 

subject to the feedback linearized state-dynamics of the system given in (70) . Given the 
relation between the linearized state and true state of the system, the corresponding gain 
matrices R and Q in (85) are chosen to minimize the appropriate control and state errors 
as 


lim z 3 (t) = 0 3 (t) = 0 nf 

/—»oo 7 J 

lim z 6 (t) = (t) = 0J ref 
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Q = diag 


'i/| A r_r.i/|* u r.i/(A*_)V 

v ! /|| Av max f , 1 /1Av^ f, 1 / (A co m J 2 y 


R = diag (l / (F„ f ,1 / (f, 1 / (7- c „_ ) 2 ) 


( 86 ) 


where r cMG max is taken to be the maximum possible imparted force and torques from the 
thrusters and MSGCMG. 

Given the use of discrete cold-gas thrusters in the system for translational control 
throughout a commanded maneuver and rotational control when the continuously acting 
momentum exchange device is unavailable, Schmitt trigger switching logic is imposed. 
Schmitt triggers have the unique advantage of reducing undesirable chattering and 
subsequent propellant waste nearby the reference state through an output-versus-input 
logic that imposes a dead zone and hysteresis to the phase space as shown in Figure 8. 



Figure 6. Schmitt Trigger Characteristics with Design Parameters Considering r BI , 

Coordinate Control Logic 


Three separate Schmitt triggers are used with the design parameters of the Schmitt 
trigger shown in Figure 8 (as demonstrated for the r BI , coordinate control logic). In the 

case of the two translational DoF Schmitt triggers, the parameters are chosen such that 


£ on( r ) = K r r db +K v V L 
£ off{ r ) = K r r db- K v V L 


(87) 
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where v, = F [mx At/2m. r db is a free parameter that is constrained by mission 
requirements. v out is chosen such that the maximum control command from the 
decoupling control law yields a value less than or equal to F max for the translational 
thruster. 

In the case of the rotational DoF Schmitt trigger when the momentum exchange 
device is unavailable, the parameters are chosen such that 


£ on K g 0 db + K (0 (O l 

£ off = KgO db ~ K co (O l 


( 88 ) 


where co L = F max LAt/2J 3 . 6 db is a free parameter that is again constrained by mission 

requirements. For both modes of operation (i.e. with or without a momentum exchange 
device), (72) through (74) can be used to determine that 

v. =v, = F / yflm (89) 

l,max 2, max max/ v k ' 


and when the thrusters are used for rotational control 

^3, max — ^max^/^3 


(90) 


When the momentum exchange device is available, the desired torque, as 
determined by the LQR control law described above, is passed directly through the 
Schmitt trigger to the decoupling control law to determine the required gimbal rate 
command to the MSGCMG. 

The three Schmitt trigger blocks output the requested control inputs along the 
7 frame. The appropriate feedback linearizing control law is then used to transform 
these control inputs into requested thrust, thruster angle and MSGCMG gimbal rate along 
the 7 frame. From these, a vector of specific actuator commands are formed such that 

=[ F a^a7 h ,a h ,ScM G ] ( 91 ) 
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Each thruster command is normalized with respect to F max and then fed with its 

corresponding commanded angle into separate Pulse Width Modulation (PWM) blocks. 
Each PWM block is then used to obtain an approximately linear duty cycle from on-off 
actuators by modulating the opening time of the solenoid valves (Wie 1998, p. 455). 
Additionally, due to the linkage between the thruster command and the thruster angle, the 
thruster firing sequence is held until the actual thruster angle is within a tolerance of the 
commanded thruster angle. Furthermore, in order to reduce over-controlling the system, 
the LQR, Schmitt trigger logic and decoupling control algorithm are run at the PWM 
bandwidth of 8.33 Hz. From each PWM, digital outputs (either zero or one) command 
the two thrusters while the corresponding angle is sent via RS-232 to the appropriate 
thruster gimbal motor. 


A6 

max 

1.8 x 10' 2 rad 

A <x> 

max 

1.8 x 10' 2 rad-s' 1 

T 

1 CMG, max 

0.668 Nm 

K r =K LQR (l,l) = K LQR (2,2) 

15.9 

K v — K lqr (1,4) = K lqr (2, 5) 

84.54 s 

K 9 = K lqr ( 3,3) 

1.39 

K a - K lqr (3,6) 

1.75 s 

r db 

10‘ 2 m 

V L 

3.05 x 10' 5 m-s' 1 

o db 

1.8 x 10' 2 rad 

co L 

1.8 x 10' 2 rad-s' 1 

£ ,J r ) 

1.61 x 10 1 m 

£ off (r) 

1.56 x 10 1 m 

£ J0) 

2.47 x 10' 2 rad 

e off 

2.37 x 10' 2 rad 

PWM min pulse width 

10' 2 s 

PWM sample time 

1.2 x 10 1 s 


Table 4. Values of the Control Parameters 
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Table 3 and Table 4 list the values of the control parameters used for the 
experimental tests reported in the following section. In particular, ||Av max || is chosen 

based typical maximum relative velocities during rendezvous scenarios while A# max is 
taken to be 1 degree and A<y max is chosen to be 1 degrees/sec which correspond to typical 

slew rate requirements for small satellites (Roser and Schedoni 1997; Lappas, Steyn, and 
Underwood 2002). The minimum opening time of the PWM was based on experimental 
results for the installed solenoid valves reported in (Lugini and Romano 2009). 

F. EXPERIMENTAL RESULTS 

The navigation and control algorithms introduced above were coded in 
MATLAB®-Simulink® and run in real time using MATLAB XPC Target™ embedded 
on the SRL’s second generation spacecraft simulator’s on-board PC-104. Two 
experimental tests are presented to demonstrate the effectiveness of the designed control 
system. The scenario presented represents a potential real-world autonomous proximity 
operation mission where a small spacecraft is tasked with performing a full 360-degree 
circle around another spacecraft for the purpose of inspection or pre-docking. These 
experimental tests validate the navigation and control approach and furthermore 
demonstrate the capability of the robotic spacecraft simulator testbed. 

1. Autonomous Proximity Maneuver Using Vectorable Thrusters and 
MSGCMG Along a Closed Circular Path 

Figure 7, Figure 8, and Figure 9 report the results of an autonomous proximity 
maneuver along a closed circular trajectory of NPS SRL’s second generation robotic 
spacecraft simulator using its vectorable thrusters and MSGCMG. The reference path for 
the center of mass of the simulator consists of 200 waypoints, taken at angular intervals 
of 1.8 deg along a circle of diameter lm with a center at the point [2.0 m, 2.0 m] in the 
ICS, which can be assumed, for instance, to be the center of mass of the target. The 
reference attitude is taken to be zero throughout the maneuver. The entire maneuver lasts 
147 s. During the first 10 s, the simulator is maintained fixed in order to allow the 
attitude Kalman filter time to converge to a solution. At 10 s into the experiment, the 


47 



solenoid valve regulating the air flow to the linear air bearings is opened and the 
simulator begins to float over the epoxy floor. At this point, the simulator begins to 
follow the closed path through autonomous control of the two thrusters and the 
MSGCMG. 

As evidenced in Figure 7a through Figure 7d, the components of the center of 
mass of the simulator as estimated by the translation linear quadratic estimator are kept 
close to the reference signals by the action of the vectorable thrusters. Specifically, the 
mean of the absolute value of the tracking error is 1.3 cm for A r B11 , with a standard 

deviation of 9.1 mm, 1.4 cm mean for Ar fl/ 2 with a standard deviation of 8.6 mm, 2.4 
mm/s mean for Av m , with a standard deviation of 1.8 mm/s and 3.0 mm/s mean for 
Av b/2 with a standard deviation of 2.7 mm/s. Furthermore, the mean of the absolute 
value of the estimated error in r BI , is 2 mm with a standard deviation of 2 mm and 4 mm 
in r BI 2 with a standard deviation of 3 mm. Likewise, Figure 7e and Figure 7f demonstrate 

the accuracy of the attitude tracking control through a comparison of the commanded and 
actual attitude and attitude rate. Specifically, the mean of the absolute value of tracking 
error for A6* 3 is 0.14 deg with a standard deviation of 0.11 deg and 0.14 deg/s for 

A<z> 5 w ith a standard deviation of 0.15 deg/s. These control accuracies are in good 
agreement with the set parameters of the Schmitt triggers and the LQR design. 

Figure 8a through Figure 8d report the command signals to the simulator’s 
thrusters along with their angular positions. The commands to the thrusters demonstrate 
that the Schmitt trigger logic successfully avoids chattering behavior and the feedback 
linearized controller is able to determine the requisite thruster angles. Figure 8e and 
Figure 8f show the gimbal position of the miniature single-gimbaled control moment 
gyro and the delivered torque. Of note, the control system is able to autonomously 
maneuver the simulator without saturating the MSGCMG. 
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Transversal CoM Position 


0.025 


Transversal CoM Velocity 



Z-axis Attitude 




Z-axis Attitude Rate 




Figure 7. Logged data versus time of an autonomous proximity maneuver of NPS 
SRL’s 3-DoF spacecraft simulator along a closed path using vectorable thrusters and 
MSGCMG. The simulator begins floating over the epoxy floor at t = 10 s. 
a) Transversal position of the center of mass of the simulator in I b) Transversal 
velocity of the center of mass of the simulator in ICS c) Longitudinal position of the 
center of mass of the simulator d) Longitudinal velocity of the center of mass of the 

simulator e) Attitude f) Attitude rate 
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Figure 8. Control actuator actions during autonomous proximity maneuver of NPS 
SRL’s 3-DoF spacecraft simulator along a closed path using vectorable thrusters and 
MSGCMG. a) Thruster 1 firing profile b) Thruster 1 position c) Thruster 2 firing profile 
d) Thruster 2 position f) MSGCMG torque profile e) MSGCMG gimbal position 

Figure 9 depicts a bird’s-eye view of the spacecraft simulator motion. Of 
particular note, the good control accuracy can be evaluated by the closeness of the actual 
ground-track line to the commanded circular trajectory and of the initial configuration of 
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the simulator to the final one. The total AV required during this experimental test was 
0.294 m/s which correspond to a total impulse of 7.65 Ns. 



1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 3 

X c (m) 

Figure 9. Bird’s-eye view of autonomous proximity maneuver of NPS SRL’s 3-DoF 
spacecraft simulator along a closed path using vectorable thrusters and MSGCMG 

2. Autonomous Proximity Maneuver Using Vectorable Thrusters Along 
a Closed Circular Path 

Figure 10, Figure 11, and Figure 12 report the results of maneuvering the 
spacecraft simulator along the same reference maneuver as in Section 6.1 but by using 
only the vectorable thrusters. This maneuver is presented to demonstrate the 
experimental validation of the STLC analytical results. As before, during the first 10 s, 
the simulator is not floating and kept stationary while the attitude Kalman filter 
converges. 

The tracking and estimation errors for this maneuver are as follows with the 
logged positions, attitudes and velocities shown in Figure 10. The mean of the absolute 
value of the tracking error is 1.4 cm for A r m ,, with a standard deviation of 8.5 mm, 1.4 
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cm mean for Ar m , with a standard deviation of 8.6 mm, 2.5 mm/s mean for Av BI , with a 
standard deviation of 1.9 mm/s and 3.1 mm/s mean for A r B[ 2 with a standard deviation of 
2.8 mm/s. The mean of the absolute value of the estimated error in r BI , is 3 mm with a 
standard deviation of 3 mm and 4 mm in r Bl 2 with a standard deviation of 5 mm. The 
mean of the absolute value of tracking error for A0, is 0.52 deg with a standard deviation 
of 0.31 deg and 0.24 deg/s for A co 3 with a standard deviation of 0.20 deg/s. These control 

accuracies are in good agreement with the set parameters of the Schmitt triggers and LQR 
design. 

Figure 11 reports the command signals to the simulator’s thrusters with the 
commands to the thrusters again demonstrating that the feedback linearized controller is 
able to determine the requisite thruster angles to take advantage of this fully minimized 
actuation system. Figure 12 depicts a bird’s-eye view of the motion of the simulator 
during this maneuver. The total AV required during this experimental test was 0.327 m/s 
which correspond to a total impulse of 8.55 Ns. 
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c) 
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Transversal CoM Position 
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Longitudinal CoM Position 



t (sec) 


Z-axis Attitude 



b) 
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Transversal CoM Velocity 



Longitudinal CoM Velocity 
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Figure 10. Logged data versus time of an autonomous proximity maneuver of NPS 
SRL’s 3-DoF spacecraft simulator along a closed path using only vectorable thrusters. 
The simulator begins floating over the epoxy floor at t = 10 s. a) Transversal position of 
the center of mass of the simulator in the I b) Transversal velocity of the center of 
mass of the simulator in ICS c) Longitudinal position of the center of mass of the 
simulator d) Longitudinal velocity of the center of mass of the simulator e) Attitude 

f) Attitude rate 
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Thruster 1 Angle 



Figure 11. Control actuator actions during autonomous proximity maneuver of NPS 
SRL’s 3-DoF spacecraft simulator along a closed path using only vectorable thrusters, a) 
Thruster 1 firing profile b) Thruster 1 position c) Thruster 2 firing profile d) Thruster 2 

position 
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Figure 12. Autonomous proximity maneuver of NPS SRL’s 3-DoF spacecraft 
simulator along a closed path using only thrusters 
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III. VECTORABLE THRUSTER CONFIGURATION FOR ROTO- 
TRANSLATIONAL CONTROL OF SPACECRAFT 

A novel configuration of two oppositely-mounted vectorable thrusters is 
presented to provide six DoF roto-translational control of a spacecraft. Figure 13 depicts 
a potential vectorable thruster solution based on the research presented by Canfield and 
Reinholz (1998). By allowing each thruster to be vectored within a hemispherical space, 
small-time local controllability of underactuated spacecraft is achieved. This 
configuration of control actuators is proposed as a possible alternative to conventional 
fixed thruster-based architectures on rapidly deployable, low-cost, and low-mass 
spacecraft systems. It may also be used for more traditional spacecraft such as NASA’s 
Crew Exploration Vehicle (CEV) and ESA’s Automated Transfer Vehicle (ATV). The 
CEV’s fixed thruster-based Reaction Control System (RCS) consists of a configuration of 
four pods of six thrusters to provide six DoF roto-translation control (Jackson and 
Gonzalez 2007). In comparison, the ATV achieves six DoF control via 32 fixed thrusters 
with 28 dedicated to the Attitude Control System (ACS) and 4 dedicated to the Orbital 
Control System (OCS) (Cavrois, Reynaud, Personne, Chavy and Strandmoe 2008). The 
potential advantages of reducing the number of thrusters required for 6-DoF roto- 
translational control include simplification of the overall spacecraft design, reduction of 
the required fuel for a desired maneuver, and a reduction in the overall size of the 
propulsion system. 
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Figure 13. Depiction of potential hemispherical vectorable thruster (Courtesy Dr. 
Steven Canfield, Tennesse Tech University) 


This chapter presents the dynamics model of this novel unconventional spacecraft 
architecture within a proximity operations environment and, following the mathematical 
methods presented in the previous chapter, the small-time local controllability of this 
underactuated system is demonstrated. By adding two additional pairs of fixed thrusters, 
the system becomes fully actuated and thus input-output linearizable. Given this input- 
output linearizability, a feedback linearized control law is derived to that can control the 
6-DoF roto-translation problem via traditional linear control methods. The reference 
spacecraft is assumed to be rigid with principal moments of inertia aligned with the body- 
fixed control axes. 


A. SPACECRAFT RELATIVE MOTION DYNAMICS 

Consider two close-orbiting spacecraft moving in the gravitational field of Earth 
or one chaser spacecraft moving with respect to a reference point. It is assumed that the 
chaser-target pair are homogeneous rigid bodies and that the chaser-target pair are 
moving in a nearly circular orbit with the final desired state of the control problem being 
a soft dock of the two vehicles. 
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Three reference frames are established: an inertial frame I with orthogonal axes 
defined by the set of unit vectors j/,, i 2 , i 3 j, a rotating satellite reference frame Jf with 

orthogonal axes defined by the set of unit vectors j /?,, li 2 , h } | and a body-fixed 

framed with orthogonal axes defined by the set of unit vectors \b v b 2 ,b 2 \. The rotating 

satellite reference frame J-f is used to define the reference coordinate system for the 
chaser-target rendezvous problem where the is collinear with the position vector of the 

target spacecraft with respect to the Earth’s center, h 2 is in the direction of the velocity 
vector of the target spacecraft aligned with the local horizontal and h, is normal to the 
orbit plane. In this manner, the angular velocity of 3~f with respect to I is 

CO/// tO qItj (92) 

where <o 0 is the target spacecraft’s mean motion. The angular velocity of 2> with respect 
to I is then given by 

(® B1 — ^HI ~ ®BH ~ (93) 


1, Translational Motion of the Chaser-Target Pair 

The 3-DoF translational motion of the rigid chaser spacecraft with respect to Jf 
are given by Vallado (2001): 


r BH ~' V BH 

VBH = f {^BH ’ V BH ) + ^ 


(94) 


wherev(t)eM 3 ,Vv = [v B // 1 ,V/ j // 2 ,V/ j // 3 ] is the translational velocity vector of the 
chaser’s body-fixed frame J 7 with respect to the moving satellite reference frame Jf 
expressed in Jf , r(t) e M 3 ,Vr = \_r BH1 ,r BH 2 ,r BH<3 ~\ is the position vector of the chaser’s 
body-fixed frame 2? with respect to the moving satellite reference frame 2> expressed 
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in F> , "F(f) e R 3 , V H F = ["Fj, H F 2 , W F 3 ] is the control vector acting on J~f , m is the 
mass of the chaser spacecraft and f (r BH , \ BH ) e R 3 represents the drift vector 


f (^BH ’ ) 


+ ^‘ (a 0 V BH,2 

—2 co Q v BH A 
—oyr 

0 'BH,3 


(95) 


2. Rotational Motion of the Chaser Spacecraft with Respect to the 
Inertial Frame 

The rotational motion of the chaser spacecraft about body-fixed axes with their 
origin located at the body’s center of mass can be described using Euler’s equations. For 
the general case, in which the control axes do not coincide with the principal axes of 
inertia and torquing devices such as independent gas thrusters are used to impart torques 
about these control axes, these become from Wie (1998, p. 341): 

<h BI = -J l a) BI J(o B1 + J l T (96) 

where <o BI (t) e R 3 , Vto BI =[o B/1 ,o B/2 ,o B/3 ] is the angular velocity vector of the body- 
fixed frame 2B with respect to an inertial frame I expressed in 2 >, 
T(l)eR 3 , VT = [ B 7j, B T 2 , B T 3 J is the control torque vector acting on 2B. The notation 
0 ) BI e R 3 ' 3 denotes the skew-symmetric matrix 


co B i 


0 


-eo B 


co E 


co, 


BI ,3 


CO 


BI, 1 


-CO 


BI, 2 


CO 


BI, 1 


0 


and J 


D 3.v3 


is the inertia matrix with respect to the JF given by 


J = 


Ji i J 


12 


J 2 1 J 22 J 23 


J 3\ J 32 


33. 


(97) 


(98) 
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3. 


Quaternion Kinematic Representation of the Chaser Spacecraft’s 
Orientation 


The orientation of a rigid spacecraft can be defined through various 
parameterizations of the special orthogonal group 50(3). In this work, the quaternion is 

considered to parameterize this orientation due to their suitability for onboard real-time 
computation and subsequent common use for spacecraft attitude determination (Wie et al. 
1989). The unit quaternion, first envisioned by Hamilton in 1843, can be separated into 
two parts, the first being a three-dimensional vector indicating the direction of the axis of 
rotation and the second being a scalar quantity which relates the angle of rotation about 
this axis of rotation such that 


q=q + q 4 =q ] i + q 1 J + q 3 k + q 4 (99) 

where q(r) = |q(r),^ 4 (r)| e R 3 xR and 

q = b sin(ri/2), i = 1,2,3 

V ’ (100) 

q A =cos(ri/2) 

where 5 represents the magnitude of the rotation around the chosen axis and 
(b\,b 2 ,b 3 )are the direction cosines locating the axis of rotation to the inertial frame by 
Euler’s Theorem (Wie and Barba 1985; Wie et al. 1989; Cristi and Burl 1993). The 
properties of the unit vectors (i, j, k)with respect to the chosen reference frame when 
multiplied are such that 0 

ij = -ji = k, jk = -kj = i, ki = -ik = j,ii = jj = kk = -1 (101) 

By considering the direction cosines expressed in (100) to be equivalent to the 

unit vectors (i, j, k), the quaternion can be used to define a sequence of rotations about a 
series of Euler angles defined by (j)J), and (// . By considering a 3-2-1 body-fixed 

rotation sequence from to I . as is typical for visualizing spacecraft orientations, a set 

of elementary quaternions can be defined from the three Euler angles similar to the three 
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sets of elementary rotation matrices (Wie 1998). In this manner, 8 X = (j),S 2 = 8,S 3 =y/ and 
from (100) and (101), the quaternion associated with this sequence of rotations becomes 


1 

1_ 


W 3 “ W 3 

9b/,2 


C,5 2 C 3 + S l C 2 S 3 

9b/,3 



1 


CjC 2 C 3 + .V, ,S’ 2 ,S’ 3 


( 102 ) 


where s f = sin(4/2), c t = cos(c> ; ./2),/ = 1,2,3. Note that the order of elementary 

quaternion multiplications is reversed from the order of elementary rotation matrices due 
to the properties of quaternion multiplication expressed in (101). 

The kinematic differential equation of the quaternion representing the orientation 
of the chaser spacecraft’s body-fixed frame 2? to I is given as 


9b; 2 (9b/,4^3*3 9b/)®B/ 

• __I 7- 

9 b /,4 


(103) 


where q' BI e M 3v3 represents the skew-symmetric matrix 


9b/ 


0 — 9b/, 3 

9b/, 3 0 

- 9b/,2 9b/,1 


9b/,2 

— 9 b/, i 

0 


(104) 


and 7 3v3 represents a 3x3 identity matrix. One can quickly show using (100) that the 
quaternion is subject to the constraint 

9b/,4 — 1 — 9b/9b/ ~ 1 — |Qb/ I” (105) 
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4. Angular Motion of the Chaser Spacecraft with Respect to the Satellite 
Reference Frame 

Considering the rotating satellite reference frame 3~f , the quaternion describing 
the orientation of 22 with respect to Jf is defined as above so that 
q, ;// (t) = |q w;/ ( t),q BH 4 (f)} e R 3 xR . From (103), the kinematic differential equation of 

4 bh is § iven b y 


q BH ~ {fiBH ,4^3x3 


C lBHA 2 q BH^BH 


(106) 


Likewise, the constraint on the components of q BH hold so that (105) becomes 


2bha i q^ 


(107) 


The 3x3 rotation matrix that brings J~f onto 22 is given by 


R-h [ c 1bha q bh ^3x3 + ^q Bif q BH 2 c[bha c Ibh 


(108) 


so that 


1 

JS"> 
_1 


1 

_1 

- 

4 

= h r b 

4 

= 

1 

1_ 


_1 

- 


R 2X 


R 


31 


^12 ^22 ^32 

i? 23 i? 33 


(109) 


From (93) and (108), the angular velocity of the body-fixed frame 22 with 
respect to 2f expressed in 22 is therefore 


^ HI + R h [O, 0, G) 0 ] 


(HO) 


where w BH (t) € M , Vt o BH — ~ co HH 1, &>bh ,2 > ^rh ,3 ] 


Taking the time derivative of (110) yields 
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(Ill) 


which, given 


where co x BH e R 3x3 


®BH ^h[0;0, (U Q ] 


B ^[0,0,fflbf =a) x BH R u [0,0,(O 0 ] 7 


represents the skew-symmetric matrix 


0 3 ^BH.2 

®BH,3 0 —®BH, 1 

~ < ^ > BH,2 ®BH,l 0 


yields 


(O 


BH 


— o> /;/ + <x> m R H 


[0,0,ruj 


( 112 ) 


(113) 


(114) 


Therefore, by substituting (96) into (114), the angular motion equations of the 
chaser spacecraft’s body-fixed frame S with respect to the rotating satellite reference 
frame J~f expressed in 2? become 


QbH 2 X^BHA^xd +< ?Bh)®BH 

1 r 

*?BH,4 = — ^Qbh^bh (115) 

d> BH T (O bi J(O bi G)qG) BB [ [-^13, -^23 ’ ~^33 ] ^ T 


5. Chaser Spacecraft Relative Motion Dynamics 

From (94) and (115), the full set of equations describing the relative motion of the 
chaser spacecraft with respect to the rotating satellite reference frame 3~f become 
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r BH V BH 

' BH = f ( l "/>« ’ V BH ) + ^/> F 

QbH = ^(^BH,4^3x3 +C lBH )®*BH (H6) 

' = _I r 

^ BH A 2 ^ ® «W 

BH ^ bI ^O^BH [^13 ’ ^23 ’ ~^33 ] J T 


with F(t)el 3 , VF = [ ~_ B F r B F 2 , B F ,J representing the control thrust vector acting on the 
body-fixed frame 2> and H R B = B R T H representing the 3x3 rotation matrix that brings 
2? onto J~f. 

Following Wie (1998, p. 407), the system of concern can be reduced to the twelve 
dimensional space through substitution of the constraint on the components of the 
quaternion given by (105) into (108) and (116) to yield 

*BH ~ V BH 

V BH = f ( *"«// ’ V IUI ) + m RIS ^ 

1/ /—-^ \ (117) 

— "2\ ’ _ I h*3 + QbH J®BH 

0) B H ^BIBl ^ 3 ) ^BH \ 3 ’ "^23 ’ "^33 ] ^ ^ 


where 


X = (i - 2 has f ) 4 v 3 + 2q BH^BH + 2 V 1 -||q i 


Qbh 


(118) 


By considering a diagonal zed inertia matrix so that (98) simplifies to 


J = 


J n 0 0 

0 J 22 0 

0 0 J 33 


(119) 


the relative motion equations described by (117) can be compactly expressed in proper 
control-affine form as 
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x = f (x) + G(x)u 


(120) 


where x(l)eR p is the 12 -dimensional state vector given by 

X — [x i ,X 2 ,X 3 ,...,X 12 ] — \^> BH p r BH,2’ } BH, 3 ’QbH .vQbH ,2’QbH ,3 ’ V BH ,1 ’ V BH ,2 ’ V BH ,3 ’ ®BH ,2 ’ ®BH ,3 ] (121) 

u(r) e U 6 is the 6 -dimensional control vector given by 


u = [mj, u 2 , w 6 f = [ B F l , b F 2 , b F 3 , %, %, %] 7 


( 122 ) 


the drift vector f(x)el p given by 


f(x) = 


— — x [Jt + a» o ( 2r;x i + 2xx )] + — I" x + ® o (I - 2x~ - 2;t~ )1 + — // [ jc + o) (2xx-2,ix)\ 
2 2 2 

~ x 6 [* 10 + ®„ ( 2 x 4 x 6 ~ lr l x <)\ + - x X x n + ® 0 (l-2x 4 -2x;)] + -/;[xii + <« 0 (277^ + 2x X 6 )] 
2 2 2 

— * 5 k + ( 2x . x e ~ 2, J X 5 )] + ~ x 4 [*„ + ( 2 1 x i + 2x s x , )] +~ T l\x n + co 0 (l - 2x - 2x )] 


3 ax + 2 co x 
0 1 0 8 


-2 co x 

0 7 




(^,3-4.) 


[x u + co a (l - 2jt“ - 2x s [*, + co a (2 rjx t + 2 xx a )] + a> o [y | (l - 2x' t - 2y ) - x u (2 r/x t + 2yy)] 
\x, 2 + ® (l - 2x a - 2x ~)] [y o + ® (2.v 4 y - 2 r/x s )] + ® \x, 2 ( 2x 4 x 6 ~ 2 t 7 x , ) - x m (l - 2x* - 2x ~)] 


(j -J ) 

\ B,1 B ,2 / 


[*n + ^0 (2t]x a +2xxJ\[ x + a> o (2x a x s - 2ijx r )] + ® o [x o ( 2i)x t + 2xx ) - x n (2x t x 6 - 2j)x a )] 

(123) 


where 77 = ±*/l-- q BH and the control matrix G(x) e R 12a6 given by 


66 



(124) 


G(x) 


J 6x6 


G, (!) 


where 6) (x) e R 6 ' 6 represents the control distribution matrix given the placement of the 

control inputs with respect to the spacecraft’s mass center. For the ideal case, where each 
thrust vector is parallel to its corresponding body-fixed axis acting through the 
spacecraft’s mass center and three independent paired thrusters provide torque about their 
corresponding body-fixed axes, Gj(x) becomes 


Gi(x) 


m ,,R B °3*3 

0 3x3 J 1 _ 


(125) 


B. SMALL-TIME LOCAL CONTROLLABILITY CONSIDERATIONS FOR 
THE PROPOSED 6-DOF SPACECRAFT DESIGN 

Small-time local controllability for the 6-DoF nonlinear equations representing 
the relative motion of a chaser spacecraft given by (120) can be determined following the 
same methods as presented in Chapter II for the 3-DoF spacecraft simulator. In order to 
properly deal with the momentum of the drift vector term f(x), the analysis of small¬ 
time local controllability is performed at the equilibrium point corresponding to f (p) = 0 
or 


P 


[0,x 2 ,0,0,0,0,0,0,0,0,0,0] 7 


(126) 


By varying the control vector u, the minimum number of control actuators for the 
system can be gained. In order to consider a general case, the thrust vector in the b x 
direction is defined as H F ] and is taken to be directed through the spacecraft’s mass 
center so that no torque is imparted on the body; the thrust vector parallel to the b 2 axis is 
defined as 11 F 2 rgt3 and is taken to be a distance L from the spacecraft’s mass center along 
the b\ axis so that both a force is imparted in the b 2 direction as well as a torque about 
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the b 3 axis; and the thrust vector parallel to the b 3 axis is defined as B F 3 mt2 and is also 
taken to be a distance L from the spacecraft’s mass center along the b x axis so that both 

a force is imparted in the b 3 direction as well as a torque about the b 2 axis. Additionally, 
the effects of pure control torques are considered. Each torque is assumed to provide a 
pure torque about its corresponding body-fixed axis so that 11 T } represents a torque 

about b x , b T 2 represents a torque about b 2 , and B T 3 represents a torque about b 3 . In this 
manner, a small-time local controllability analysis is conducted on variations of the 
control vector u(t) e U 6 , Vu = [ B Ei, B F 2lrot3 , B F 3lrot2 , B T X , B T 2 , B r 3 ] such that the 

corresponding vectors of the associated control distribution matrix G x (x) are: 


B F x —> g(x) = \jrF x R x ,,m 'R x2 ,m 'R u , 0,0,0] 

(127) 

% rof3 -> g( x ) = [in'R 2X ,m-'R 22 ,m-'R 23 ,0,0,LJ 3 'J 

(128) 

BF vroti -> g ( x ) = [ m l R 3X , m'R 32 , m'R 33 ,0, LJ 2 , o] 

(129) 

B r 1 ->g(x) = [o,o,o,zj 1 _1 ,o,o] r 

(130) 

B r 2 ^ g (x) = [o,o,o,o,L/; 1 ,o] r 

(131) 

^ g (x) = [o, o, o, o, o, l/ 3 1 ] r 

(132) 


The results of this investigation, which are summarized in Table 5, demonstrate 
that the minimum number of actuators for small-time local controllability of the relative 
motion spacecraft dynamics defined by (117) include one thrust vector acting through the 
spacecraft’s center of mass and two body-fixed control torques. This combination can be 
achieved by considering a total of six control actuators with two fixed thrusters mounted 
on opposite faces of the spacecraft with their thrust vectors acting in an opposing manner 

along the b x axis through the spacecraft’s center of mass and two sets of fixed paired 
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thrusters acting together to produce pure torque about two of the body-fixed axes. For 
instance, by considering the control vector u(t) e U 3 , Vu = [ S F P B T 2 , B T 3 j and the 

corresponding control distribution matrix G 1 (x) e R 6r3 


G 1 (x) = [g 1 ,g 2 ,g 3 ] 


m 1 R n 0 0 

m l R n 0 0 

m l R u 0 0 

0 0 0 

0 ZJ 22 0 

0 0 ZJ33 1 


(133) 


the Lie Algebra evaluates to: 


-Z’( A) = span < 


gpg 25 g 3 ’[ f ’gl]’[ f ’g 2 ] 5 [ f J g3]>[ f ’[ f ’g 1 ]]’[ f ’[f,g 2 ]], 

[gp[ f ’g 2 ]]{gp[ f ’g3]]’[g2’[ f ’g3]][ f ’[ f ’[ f ’gl]]] 


(134) 


However, by considering the vectorable thruster depicted in Figure 13, only two 
control actuators are required to achieve small-time local controllability of the relative 
motion spacecraft dynamics given by (117). This can be achieved by placing each 
hemispherically vectorable thruster a distance L from the spacecraft’s center of mass 

along the b x axis. For this configuration, the control vector 

u(t) g U 3 ,Vu = [ b F p B F 2/rot3 , B F 3/wt2 j is considered with the corresponding control 
distribution matrix G, (x) e R 6 ' 3 


G 1 (x) = [g 1 ,g 2 ,g 3 ] = 


m 

X 

m R 2l 

m R 3l 

m 

%2 

m l R 22 

m l R 32 

m 

~% 3 

m l R 23 

m 1 /? 33 


0 

0 

0 


0 

0 

LJ-] 


0 

ZJ3-3 

0 


(135) 


yielding a Lie Algebra for this minimally configured system of 
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£ (A) = span 


gj , g 2 , g 3 , [f, g, ] , [f, g 2 ] , [f, g 3 ] , [f, [f, g, ]] , [f, [f, g 2 ]] , 

[gi»[ f »g 3 ]]»[g 2 »[ f . S 3 ]]. [ f . [gi»[ f »g 3 ]]] * [ f > [g 2 . [ f . g 3 ]]_ 


(136) 


Control Inputs 


dim(£(A)) 

Controllability 

B 7-> B T-> B T-> 

F x or E 2/ror3 or F yrot2 


3 

Inaccessible 



6 

Inaccessible 

{‘F,,T : ) 


8 

Inaccessible 

("r.T, ),( b F„T 3 ),(%„, 2 ,T 2 ).( , f 2 „„, 
( B F s ,„ 2 ,T s )or(’F 2lm ,,T 2 ) 

,T,), 

9 

Inaccessible 

( BF 2,rot3 . "Xro.l ) > ( ^ 2 /™3 - T X ) °> ' ( ^3/r 0 ,2 

X ) 

10 

Inaccessible 

( B F l , %,%),( b F { , % , %) or ( b F { , b T 2 , 

%) 

12 

STLC 

B 77 B t? B T? 

F \ ’ r 2/rot3 ’ r 3/rot2 


12 

STLC 


Table 5. STLC Analysis for the 6-DoF Relative Motion Spacecraft Dynamics 


C. PROPOSED FEEDBACK LINEARIZABLE SPACECRAFT SYSTEM 

DESIGN WITH MINIMUM NUMBER OF CONTROL ACTUATORS 

In the previous section, a full analysis of the small-time controllability for a 
spacecraft with variations on the control inputs was conducted. This analysis yielded two 
satisfactory variations in control inputs for small-time local controllability: a set of six 
control actuators consisting of two opposing fixed thrusters acting through the 
spacecraft’s center aligned with one of the three body-fixed axes and two sets of paired 
thrusters acting about two of the three body-fixed axes or a set of two control actuators 
consisting of oppositely mounted hemispherically vectorable thrusters. From this 
analysis, a unique spacecraft control actuator configuration is proposed that leverages this 
study to provide a feedback linearizable system as might be necessary in a path or sensor 
constrained environment. As stated in Section II, in order for a system to be state- 
feedback linearizable, it must be square so that the number of actuators is equal to the 
number of observed states. From this, the number of observed states and number of 
control inputs are assumed to be six and a novel minimally actuated system is proposed 
consisting of a total of four thrusters. Figure 14 graphically depicts the forces acting 
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upon the body-fixed frame 2? through two hemispherically vectorable oppositely 
mounted thrust vectors F fl e 1R 3 and F /; e 1R 3 and a set of paired thrusters acting about the 

/?, axis to impart the torque fei. 



k 


Figure 14. 6-DoF Minimally Actuated Feedback Linearizable Spacecraft 

Configuration 

In order to properly resolve F a and F /; into 2 ?, the vectorable thruster frame 
Tfis established with orthogonal axes defined by the set of unit vectors {v (l ,v j2 ,v )3 } 
where i = a,b corresponding to either F fl or F b . Both axes v a2 and v a2 are taken to be 
collinear with the b 2 axis. Additionally, the v a3 axis is taken to be aligned with the b x 

axis with v fl3 = b t while the v h2 axis is taken to be aligned with b x with v /)3 = —b x . The 
coordinates of the imparted force vector, as depicted in Figure 14, is defined in the 
thruster frame through two angles a i and fi , referred to as the polar and azimuthal 

angles, respectively. The azimuthal angle a,, is the angle between the axis of the 
imparted thrust vector and the v (3 unit vector and is limited to the range [0, zr/2]. The 
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polar angle is defined as the angle between the v n and the projection of the thrust vector 
F ( . onto the {v n ,v i2 } plane. In this manner, 


' F „ = ~K [ sin a a cos P a > sin a a sin Pa > cos cc a f 


and 


' F b =-F b [sin a h cos J3 b , sin a b sin J3 h , cos a h ] 7 


The 3x3 rotation matrices that bring 7f onto 2> are then 



0 0 1 

0 1 0 

-10 0 


and 



0 0-1 
0 1 0 

1 0 0 


(137) 


(138) 


(139) 


(140) 


In this manner, the hemispherically vectorable thrust vectors can be expressed in 
the body-fixed frame 2? as 

B = ~F a [cos a a , sin a a sin fi a , - sin a a cos J3 a ] r (141) 

and 

B F b = -F b [-cosa b ,sina b sin j3 b ,sina b cos fi b ] T (142) 

yielding the thrust vector given in (117) by F (t) e R 3 , VF = /f F ; + 11 F /; . 

The torque vector in (117) given by Te M 3 can be resolved as 
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(143) 


T = 


T c 

~F a sin cc a cos p a - F h sin a b cos p b 
~F a sin cc a sin p a + F h sin a h sin p b 


Defining the control inputs as such yields the control vector 

u (?) e U 6 , Vu = [ B F l , B F 2lrot3 , B F 3lrot2 , >t T ] , b T 2 , and corresponding control distribution 

matrix G 1 (x) e R 6 ' 6 


m 

X 

m a 2I 

m R, 3 

0 

m Y R 3l 

m R 2l 

m 

X 

m R 22 

m R 32 

0 

m R 32 

m R 22 

m 

'% 3 

m l R 23 

m l R 33 

0 

m l R 33 

m l R 23 


0 

0 

0 

dJn 

0 

0 


0 

0 

dJ 2 2 1 

0 

~dJ~ 22 

0 


0 

-dJ 33 

0 

0 

0 

dJ 33 


The components ^represent the entries in the corresponding i-th row and /-1h 
column of H R B as defined by (109) and the variable d is defined similar to what is 
proposed for the 3-DoF spacecraft simulator presented in Chapter II whereby the control 
variables are defined given a desired control input to the system with respect to b x . Thus 
if 


B C, < 0 -> 


B C, > 0 -» 


d = —L 

* F i = -F a cos a tt , B F 2/rot3 = -F a sin sin p a , B F Vrot2 = F a sin a tt cos p a 
% = T c , b T 2 = -F b sin a h cos p h , % = -F b sin a b sin p b 

d = L 

b F x = F b cos a b , B F 2lrol3 = -F b sin a h sin p b , B F 3lrot2 = -F b sin a h cos p b 
X = -T c , X = F a sin a a cos p a , B T 3 = -F a sin a a sin p a 


(145) 


1. Smooth Feedback Control via State Feedback Linearization 

The methodology presented in (62) through (65) can be used to demonstrate that 
the nonlinear relative motion dynamics given by (117) are state feedback linearizable. 
With the nonlinear relative motion dynamics expressed in control-affine form given by 
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(120) with G, (x) defined by (144), the matrix £(x) as defined by (63) is invertible and 

the relative degree of the system at the equilibrium point p is 

(r 1 ,r 2 ,r 3 ,r 4 ,r 5 ,r 6 ) = (2,2,2,2,2,2). Therefore, the total relative degree of the system is 

twelve. Given that the total relative degree of the system is equal to the number of states, 
the nonlinear system can be fully linearized by state feedback to yield smooth feedback 
control. If we choose 

h ( x ) = [ h i( x )’ h 2 ( x )T =[[xi,x 2 ,x 3 ],[x 4 ,x 5 ,X6]] r =[r BH ,q BH ] T (146) 
the state transformation can be chosen as 


z = [hj ( x ),h 2 (x), L { h, (x),L t h 2 (x)] r 




(147) 


That is 


r BH ~ Z 1 


*\.BH ~ Z 2 


Z 3 


(148) 


01 BH ~ Q Z 4 

o) BI =Q\ + b R h [0,0,co 0 ] T 

where z l eR 3 ,z 2 eR 3 ,z 3 eR 3 ,z 4 e R 3 are new state variables and Q e M 3x3 is given by 


Q = 


2 ( — V^ II ^3x3 + QbH 


(149) 


and the system is transformed into 
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(150) 



z 3 =f(z 1 ,z 3 ) + H i? 5 F 

J l S (Q \z 4 + b R h [0,0, co, f ) j(Q \ z 4 + b R h [0,0, co, f ) 
[5(G~ 1 z 4 )X[0,0,q ) f +7 _1 T 


z 4 = QQ~\+Q 


+ 


where 5 (•) represents the skew-symmetric matrix of • . 

The relative motion dynamics given by (117) can now be transformed via (150) 
and the state feedback control law 

[F,Tf =£(x) -1 (v-b) (151) 


into the linear system 


^6x6 

16x6 

z + 

^6x6 

_^6x6 

0 6 * 6 J 


1 


L 6x6 J 


(152) 


where £"(x) e R 6x6 = L c L f h(x),b e R 6 = Z^h(x) and v(t) eR 6 ,Vv = [v p v 2 ] r . 

With the system now in linear form, traditional linear control techniques such as 
the linear quadratic regulator as presented in Chapter II can be used to stabilize it to a 
desired state. 


2. Determination of the Commanded Azimuthal and Polar Thruster 
Angles 

Using the decoupling law as expressed in (151), the commanded azimuthal and 
polar angles for the hemispherically vectorable thrusters and the required torque for the 
paired thrusters can be determined. For example, if B U { > 0 with the azimuthal angle 
0 < a h < njl and the polar angle 0</3 b <2n, the angles to slew thruster b can be 
resolved as 
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a h = tan 


l B p2 + B f 2 
1 2/rot3 ^ 1 Vrot2 


(153) 


A = atan 2 (V 2/ro( 3 ,V 3/ro(2 ) 


(154) 


and the required thrust as 


5 77 B t?2 . B t7>2 . B 77 2 

A=v A + A/ rof 3 + A / rot 2 


(155) 


The angles to slew thruster a can be resolved as 


(156) 


A= atan2('r„*r 2 ) 


(157) 


and the required thrust as 


5 77 > / Brji2 . Brj-i 2 

A = V A + A 


(158) 


Likewise, for the case where B C/j < 0 with the azimuthal angle 0 < <7i/2 and the polar 

angle 0 < J3 b < In , the angles to slew thruster a can be resolved as 


a„ = tan 


— B F 

-1 _f_l_ 

/ 5 77 2 1 ^ 7^2 

V V r 2/rot3 " r r 3/rot2 


(159) 


A =atan2( B F 2/ra;3 ,- B F 3/ro(2 ) 


(160) 


and the required thrust as 


5 77 B y-i2 . 5 7 t>2 . 5 772 

A=V A + A/raf3 + A/n*2 


(161) 


The angles to slew thruster b can be resolved as 
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(162) 


and the required thrust as 


n 



A=atan2 (%,%) 


B 77 > / Brj! 2 i Brji 2 

= V 7 2 + 7 3 


(163) 


(164) 
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IV. QUATERNION FEEDBACK ATTITUDE REGULATION OF 
UNDERACTUATED RIGID SPACECRAFT 

Methods of providing three-axis stabilization via two body-fixed control torques 
are introduced in this chapter. The current state-of-the-art in three-axis stabilization for a 
rigid spacecraft with arbitrary inertia via two body-fixed control torques is limited to rate 
stabilization followed by three sequential Euler angle rotations or time-varying attitude 
control based on a variety of switching logics. Byrnes and Isidori (1991) demonstrated 
that for the full underactuated rigid spacecraft system of angular motion equations, a 
smooth, time-invariant state feedback control can be designed to locally stabili z e the 
spacecraft with a revolute motion about the unactuated axis. Without violating their 
assertion, this chapter introduces a novel, time-invariant smooth state feedback controller 
based on generalized inverse methods and quaternion kinematics to obtain simultaneous 
attitude stabilization to an arbitrarily small region of the zero error state with two 
bounded body-fixed control torques for an asymmetric spacecraft. This smooth state 
feedback controller will be shown to be most naturally suited to either rest to rest 
reorientation maneuvers or detumbling and reorientation maneuvers. It is not applicable 
for attitude mainentenance for attitude tracking or in the presence of large disturbance 
torques as a perturbation is imparted about the unactuated axis through the interactions 
between the two actuated axes to stabilize the attitude to a desired orientation. This 
chapter provides two key contributions: First, the method of generalized inverse control 
is extended from a study of only the dynamics as considered by Bajodah (2009b) to the 
full system of dynamic and kinematic equations. Second, two separate null-control 
vector designs are derived for the proposed quaternion feedback regulator which 
capitalize on the generalized inverse control methodology. The first design is Lyapunov 
function based and is conjectured to yield global stability but results in significant control 
effort noise near an inherent singularity in the null-control vector. The second design is 
based on perturbed feedback linearization and provides local stability while keeping the 
control effort to low levels. In addition to formal proofs for each control design, 
simulation results are provided to demonstrate the capabilities of these novel controllers. 
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A. UNDERACTUATED RIGID SPACECRAFT REORIENTATION 
THROUGH QUATERNION FEEDBACK 

In this section, the general case of an underactuated rigid spacecraft being rotated 
by two body-fixed control torques is considered. In order to simplify the problem, an 
ideal control torque is assumed. The rotational motion and quaternion kinematic 
representation of the orientation of the rigid body remain as developed in Section III for 
the general case of the body-fixed frame 22 with respect to an inertial frame I . In 
summary, these are: 


cb B/ =-J x co x bi J(3) bi + / 1 T 

Qb/ = ^(^B/,4^3.t3 (165) 

' = _I 7 

2 BIA 2 ^ BI ^ BI 

1. Attitude Error Parameterization via Quaternions 

The 3x3 rotation matrix that brings I onto 22 is given by 

^1 — (^B/,4 _ Qb/Qb/ ) ^3x3 — ^2bia2bi (166) 

Similarly, the desired attitude of the spacecraft can be described by a desired, 
body-fixed frame 22 d with respect to I . In this manner, the attitude is parameterized by 

considering a desired unit quaternion q B ^(f) = jq B fl (I4 (f)} e M 3 .vM where the 

constraint of (105) holds. Therefore, the 3x3 rotation matrix that brings I onto 22 d can 
be defined as 


R/ - 2b,,!A QBjIQBjI )^3x3 + ^QbjIQBjI ^2 B d l , 


4 2b a 


(167) 


and from (108) and (167), the 3x3 rotation matrix that brings 22 d onto 22 is therefore 


RB d ^1 B,BA 0 B,, bO B,,.B ) ^3x3 + ^1 Hi /;*! B ( B ^2 B,Ba2 B,B (168) 
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where the quaternion error vector q B jB (t) = jq B /B /B4 (*)} e K 3 .rM is defined by 
Behai et al. (2002) as 


q B d B QBjIA^BI 

T 

QbjB — QbiaQbjIA +q B /q Brf / 


(169) 


The error quaternion kinematic differential equation then becomes 


q 


B d B 


QBjB 


2 (^B d BA^3x3 QBjB^^BI 

1 r 

2 q B d B & Bl 


(170) 


with (o BI = (o B B considering a rest-to-rest maneuver where ( o BjI = 0 3xl . For simplicity, 

the components of q B B will be referred to as q /; / , = {q t ,q 2 ,q 3 ,q 4 } for the remainder of 

this chapter. The full rigid spacecraft system of angular motion equations considering the 
error quaternion kinematics can be expressed as 


q B d B 2 i^B d BA^3x3 + QB d B ) 01 BI 

■ = _I r 

C lBjB 2 QbjB^BI 

c o BI — —J d) B jJ{0 + J T 


(171) 


2. Underactuated Rigid Spacecraft with Two Control Torques 

In considering an underactuated rigid spacecraft with only two independent 
control torques, the unactuated axis can be taken to be without loss of generality and 
(96) can be rewritten as 

cb B/ — F BI ^(o B i t (172) 

where x(t) e ]R 3 ,Vt = 7 _1 T = [0,u] r is the scaled control vector with 
u(t) e U 2 , Vu = [w 2 ,u 3 ] and F(to) e M 3a3 denotes the drift matrix 
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J ® bi J 


( 173 ) 


In this manner, the full underactuated rigid spacecraft system of angular motion 
equations become 



_ 1 / 

X \ 

^ B„B 

- 2 \^B d BA^3x3 



1 r 


c Ib l ,b 

2 ^ B d B®*BI 


®BI 

fc, 

II 

+ T 


(174) can be written in control-affine form as 

x = /(x) + G(x)u = F(x) & bi +Gu 


(174) 


(175) 


with state vector x(t) e M 6 vIR, Vx = [q l ,q 2 ,q 3 ,q 4 ,a> l ,a) 2 ,a> 3 f , control vector 
u(t) e M 2 ,Vu = [Mj,M 2 ] r , drift matrix F(x)eR 7x3 given by 


F(x) 


2 (^4^3x3 + C JBjB ) 

1 r 
2 ^ 

—J x co' BI J 


(176) 


and control matrix G e M 7x2 given by 


G = 


5x2 


1 2x2 


(177) 


Furthermore, (175) can be split into two coupled differential subsystems of 
concern; the first coupled subsystem combines the unactuated dynamics with the 
unactuated quaternion kinematic component, while the second considers the remaining 
two quaternion kinematic components. Both of these systems are indirectly controlled 
through the action of the angular velocity components about the actuated body-fixed 
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axes. The third subsystem considers the actuated dynamics. In fact, although global 
asymptotic stability cannot be achieved without unbounded control, the full system of 
underactuated rigid spacecraft equations of angular motion can be stabilizable to an 
arbitrarily small region about the equilibrium of the system via time-invariant smooth 
state feedback control by considering a damping coefficient to ensure the generalized 
inverse of the controls coefficient in the proposed control law does not go unbounded at 
times during the maneuver. This concept will be discussed in detail in subsequent 
sections. 

In order to compactly view the three subsystems of concern, it is useful to first 
partitioning (173) as done by Bajodah (2009b) where 


F((o b/ ) = 


| ( W B/ ) 
^21 ( W B/ ) 


^12 (®BJ ) 
^22 ( W B/ ) 


(178) 


with F n (<d b/ ) e R U1 , F n (co B/ ) e R U2 , F 2l (to B/ ) e R 2cl , F 22 (co B/ ) e R 2 ' 1 ' 2 , and then 
consider the drift matrix associated with the vector differential equations of (170) as 



>nl 

(‘W* ) 

1 F n \ 


1 

1 

to J 


1 ^22 1 

(^b.b) 

1 


(179) 


where F n (q B(B ) e M lvl , F 12 (q BjB ) e R lx2 ,F 2l (q B;B ) e M 2vl , F 22 (q BjB ) e M 2v2 . Now, by 

taking the unactuated state vector to be x u (f)eR 2 ,Vx B = [q u , co u ] where q u = q { and 

co u = co x , and the actuated angular velocity vector to be co a (t) e M 2 , Vco u = [ty 2 ,<z> 3 ] r , the 
unactuated system becomes 


x 


u 


0 

0 


F n! 


x„ + 

1 



(to B/ ) J 

_F,2 

(®B/ )_ 


(180) 
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By considering the two remaining components of the error quaternion to be 
pseudo-actuated, the pseudo-actuated error quaternion vector can be taken to be 

q p (t) e M 2 , Vq p = [q 2 , q 3 ] 7 and the actuated error quaternion subsystem then becomes 

q P =[o 2 u ^2i(q^)] x „+^22(q^)o> fl (isi) 

The third and final coupled subsystem which considers the actuated dynamics and 
explicitly contains the control vector u is therefore 

<=[0 2 ,, ^2i(« B/ )]x„+F 22 (ra B/ )(o a +u (182) 

B. LINEAR PARAMETERIZATION OF THE UNDERACTUATED RIGID 
SPACECRAFT SUBSYSTEM 

As mentioned previously, the unactuated subsystem given in (180) can be seen to 
be indirectly affected by the control vector u through w fl . By defining the scalar 

function ^(x H ): IR 2 —» R as 

= ( 183 ) 

where a> 0 represents a scalar on q u and ^(x H ) is a continuous twice differentiable 
function satisfying 

0( X B ) = °^® B =~aq u (184) 

the unactuated subsystem (180) can now be transformed into a stable, linear, second- 
order dynamic system similar to what is proposed in Bajodah (2009a and 2009b) where 

) + 2 ^( x « ) + r 2 </>{\) = 0 (185) 

with the coefficient y chosen such that (185) is stable, i.e., the characteristic equation 

s 1 +2y s + y 2 =0 (186) 
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has strictly negative-real part roots thus implying that y > 0. 

The first time derivative of ^(x H ) along solution trajectories of the underactuated 
rigid spacecraft angular motion equations given by (175) considering the control vector 
u, the control matrix G given by (177) andf (x) = F[x)( 0 BI is: 

- x - - (f (x) + Gu) - L f </>(x u ) + LJ ^)u (187) 

where L G (f)(x u ] evaluates to the 2x1 null vector. Likewise, the second time derivative of 
<!>{x u ) along the solution trajectories of the underactuated rigid spacecraft angular motion 
equations: 

h \ ^-) * - SL ^- } (f (*)+C M ) = (X„)+44■<*(*,)» (188) 

where L f ^(x i( ), L^(x i( ) andL c L f ^(x i( ) represent the Lie derivatives, or directional 
derivatives of (f>{x u ) along the direction of the vector fields defined byf(x) and 

G(Slotine 1991). With (f> and <j> given by (187) and (188), (185) can then be expressed 
as the point-wise linear form 

a r (x)u =b(x) (189) 

where the controls coefficient a(x) e K 2aI is given by 

a(x) = [L G L f ^(x„)] r (190) 

and the scalar controls load b( x ) is given by 

b(x) = -Lt<f>(x u )-2yL [ <l>(x u )-r 2 0(x u ) (191) 

The following definitions, propositions and theorem include only minor changes 

from Bajodah (2009b) to further expand them to include the full state x containing both 
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the dynamics and kinematics of the system. The proofs of these flow directly from those 
presented in Bajodah (2009b) due to the cascading nature of the controls. 

Definition 1: The desired linear unactuated dynamics given by (185) is said to 
be realizable by the rigid spacecraft angular motion equations (174) at specific values of 
x(r) if there exists a control vector u {t') that solves (189) for these values of x(r). If this 

is true for all x(f) such that x u ( t ) ^ 0 2xl , then the linear unactuated dynamics given by 
(185) is said to be globally realizable by (174). 

Proposition 1: By letting (190) be the controls coefficient relative to (183) along 
f (x) for the linear dynamics given by (185) that is globally realizable by (174), then 

a ( x ) = 0 2 *i O<f>{\d) = 0 2xl (192) 

Proof: For a vector u to exist that solves (189) at specific values of x, b(x) 
must be in the range space of a r (x) at that value of x . This is possible for any value of 
Z?(x) only if (189) is consistent. That is to say that not all elements of a T (x) vanish at 
that value of x. Therefore, the existence of a vector x* ^ 0 such that a(x*) = 0 2vl 

implies that the desired linear dynamics of (185) are not realizable at x* , which in turn 
violates global realizability of these desired linear dynamics and thus proves sufficiency. 
Necessity flows from the fact that the elements of the drift vector f (x) are multivariable 

polynomials with only quadratic elements in the components of x and that ^(x i( ) has a 
bounded x u gradient, so that a(0 6tl ) = 0 2vl . 

Definition 2: The zero actuated full-state Jacobian of the controls 
coefficient\/ 0 (x) is the square matrix that results by partially differentiating the controls 

coefficient with respect to the actuated angular rates ( evaluated at 0 ) n = 0 2cl . 


M*) 


5a r (x) 

da a 


l ®o=# 2x1 


(193) 


86 



Proposition 2: The unactuated linear dynamics given by (185) is globally 
realizable by the rigid spacecraft angular motion equations given in (174) if and only if 


det[y o (x)]*0 


(194) 


Proof: The partial derivative of (189) with respect to and evaluating the 
resulting equation at = 0 yields 


^>( x )» = 


db(x) 


dm. 


JtO =0 ? 


(195) 


Invertibility of the zero actuated state Jacobian of the controls coefficient implies 
that a control law ucan be derived for global realizability of the desired linear dynamics 
of (185) as 


u = 



db[x) 


da )„ 


Jco =0 ? 


(196) 


which proves necessity. Consider the non-linear time-varying system given by the 
equations 

z 0 =a(z 0 ) (197) 


where z 0 = [z„,z fl ] 7 , Vz [( eM,Vz a el M . From Proposition 1, global realizability of the 
desired linear dynamics of (185) implies that the origin z n =0 2xl is the unique 
equilibrium point of (197) at z u = 0 . Furthermore, since a(z 0 ) is a smooth vector field, 
it follows from Milnor’s theorem (Bajodah 2006b) that it is also a global diffeomorphism 
on M 6x1 , i.e. that it has continuous partial derivatives and an invertible zero actuated state 
Jacobian and thus sufficiency is obtained. 

Theorem 1: Given a non-singular zero actuated state Jacobian ( x ) of the 
controls coefficient a(x) = 0 2xl along f (x) = F (x)co for cdl \(t) e R 6 , the infinite set of 
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control laws that globally realize the unactuated dynamics by the underactuated rigid 
spacecraft angular motion equations given by (174) is given by 

u = u + P(x)y (198) 

where u e M 2 is 

u = a + (x)b(x) (199) 


and a + (x) e IR 2xI represents the Moore-Penrose or generalized inverse of the controls 
coefficient a(x) e R 2 ' 1 so that 


■w= 


a 


M 


\m 

0 2 ,, 


|a( x )| 56 0 
|a(x)| = ° 


( 200 ) 


and P(x) e R ’ 1 ’ represents the null-space projection matrix of a (if) such that 


p ( x ) = 4t2+ a+ ( x K( x ) ( 201 ) 

and y e M 2 ' 1 is an arbitrarily selected null-control vector. 

Proof: Satisfaction of condition (194) implies that the desired linear dynamics of 
(185) are globally realizable by the underactuated rigid spacecraft system of angular 
motion equations. From Proposition 1, this global realizability further implies that 
a(x) ^ 0 2tl at which infinite number of solutions for the point-wise linear relation (189) 

exist. Multiplying both sides of (198) by a T (x) yields 


a r (x)u = a T (x)a + (x)6(x) + a r (x)P(x)y 

= ^( x ) 

and thus (189) is recovered. Therefore, the control vector u given by (198) linearly 
parameterizes all solutions of (189) by the null-control vector y. □ 
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In this manner, the infinite set of closed-loop full underactuated rigid spacecraft 
system of angular motion equations in partitioned form become 


x = 


0 F U (<U») 

.0 


x.. + 


-^12 ) 
F l2 (or BI ) 


0) 


®2a1 F l\ ) X w + F 22 ) 01 a 


d 4 = 


0 ~b- 


x „+q P o> 


p a 


= [°2.vl F 21 (®B/ )] X „ + F 22 +« 


(203) 

(204) 

(205) 

(206) 


(198) consists of two parts. The first part, given by u, is termed the particular solution 
and acts specifically on the range-space of the generalized inverse of the controls 
coefficient a + (x). The second part, given by P(x)y , is termed the auxiliary solution and 
resides in the orthogonal complement subspace, or the null-space of the controls 
coefficient in M 3 with the null-control vector y being projected onto this space by means 
of the projection matrix P(x). 

As discussed in Bajodah (2009b), the null-control vector yeM 2xl is not fully 
arbitrary and should be chosen to yield stability of the closed loop system of equations 
given in (203) through (206). Two such null-control vectors will be presented in the 
following sections, the first is Lyapunov function based to yield global stability for the 
system given established bounds on the control and the second is perturbed feedback- 
linearizing to yield local stability and guaranteed boundedness. The development of 
these null-control vectors and their subsequent proofs of stability provide a key 
contribution to the body of knowledge with respect to underactuated rigid body attitude 
control. 
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C. DEVELOPMENT OF THE QUATERNION FEEDBACK REGULATOR 

AND STABILITY ANALYSIS 

In this section, the stability of the closed-loop system of equations given by (203) 
through (206) is addressed for two derived null-control vectors. From the definition of 
the generalized inverse of a(x) given by (200), 

lim a + (x) = co 2xl (207) 

a (x)->0 2l i 

which implies a singularity on the closed-loop stability of the system. Similarly from 
Proposition 1, if the linear unactuated dynamics (185) are globally realizable by the 
underactuated rigid spacecraft system of angular motion equations given in (175), then 

lim a (x) = ° 2d (208) 

In order to properly bound the control input to allow for the control law 
derivations that follow, the damped generalized inverse and its corresponding damped 
controls coefficient null-projector are used as presented in Bajodah (2009a). The damped 
controls coefficient generalized inverse is formulated by considering an arbitrarily small 
damping coefficient /?, to provide a bound on the generalized inverse as the squared 
Euclidean norm of the controls coefficient tends to zero. Thus, the damped generalized 
inverse of the controls coefficient a + d (x) e R 2 tl is 

<( x ) = 

where scalar /i, is a positive damping coefficient. Bounding the generalized inverse of 
a(x) given by (200) in this manner smoothes the infinite set of control laws presented in 
Theorem 1. Furthermore, 
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( 210 ) 


and 



lim |k(x)|| = 0 (211) 


and a^(x) pointwise converges to a + (x) as /?, vanishes (Bajodah 2009a). 

Likewise, the damped null-projection matrix of a(x)is modified from (209) such 

that 


^( X ) = / 2,2- a rf( X K( X ) 


(209) and (212) imply that 


^W= 


hpiim, 


l 2x2 


l ( x )|' 




A 2 


( 212 ) 


(213) 


and from (208), during the steady-state phase of response the damped null-projection 
matrix of a(x)becomes 


lim A( x ) = 4v2 ( 214 ) 

(/){ x u )-»0 7 

such that the auxiliary part of the infinite set of control laws expressed in (198) converge 
to the null-control vector y. Construction of these control laws with (209) and (212) 
yields the damped control vector u defined as 

u,/=» d +A( x )y (215) 


where u d e M M is given by 
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( 216 ) 



D. LYAPUNOV FUNCTION-BASED QUATERNION FEEDBACK 
REGULATOR 

Let the null-control vector y be chosen as 

y = ( x )[-«lAi -™ T a*d -kq T BiB <o]~ d <o a (221) 

where k,d are positive scalar gains and the damped actuated dynamics 
projector i\ + d (x)el 2 is defined as 


<( x ) 


. A 


°>IA( X )°><AA 

<A( x )«> fl <A 


( 222 ) 


where /? 2 is a positive damping coefficient on i) ( x). In this manner, substitution of (221) 
into (212) yields the class of control laws 

= «/ + A ( x ) {< ( x ) [-«! Ai (®«) - < F 22 (to B/ ) } ( 223 ) 
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which yield the closed-loop dynamical actuated subsystem 


®a ^21 ) m u + P 22 ( W i/ ) W n + U «/ 

-olu-kq] 


+p d ( x ) 


<( x ) 


T 

1 


-dto, 


( 224 ) 


Conjecture 1. Let t/ie function 0(x i( ) Le globally twice continuously 
differentiable and satisfy the condition given by (184), t/ze controls coefficient a(x) gzvezz 
by (190) and controls load b(x)given by (191) relative to ^(x H ) of the linear unactuated 
second-order dynamics given by (185) along the drift vector f (x) = P(x)co B/ of the 
underactuated rigid spacecraft system of angular motion equations given by (175). 
Furthermore, let the damped generalized inverse of the controls coefficient n) (x) be 

given by (209), the damped null-projector of the controls coefficient P, (x) be given by 
(212), and the null-control vector y be given by (221) with damped actuated dynamics 
projector vf d (x) be given by (222). If the zero-actuated state Jacobian J t] (x) given by 
(193) is non-singular, then the positive scalar feedback gains k,d e M can be chosen to 
globally stabilize the closed-loop underactuated rigid spacecraft system of angular 
motion equations given by (217) through (220) by bounded, smooth, time-invariant state 
feedback to an arbitrarily small vicinity of the equilibrium of the system defined by the 
damping coefficients fj and fi 2 . 

Proof. Let the scalar function e(x) be defined as 

e( x ) = ^(x„) + ^(x„) (225) 

where ;/is as defined for (186). Given (189) and (225) , e(x) therefore evaluates to 

e(x) = ^( x „) + ^( x „) = -7e( x )-L( x ) + a r ( x )u (226) 

Let the Lyapunov control function be defined as 
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( 227 ) 


V -~ e2 (x) + -k V,CD a +qf +q 2 + q] +(q A ~l) 

= |e 2 (x)+|r 1 (o> a +2(i-^ 4 ) 

Note that V is positive definite and asymptotically unbounded in e(x) and« fl . 

Furthermore, V can be partitioned into two positive definite functions such that 
V = V, + V 2 where 



(228) 

k~ l G> T a (O a +2(l-q 4 ) 

(229) 


where k 1 e R is a positive scalar constant gain. 

The corresponding time-derivative of V l is given by 

V x =e , (x)e , (x) = -e(x)7>(x)-7e 2 (x) + £ , (x)a r (x)u (230) 


which by is reduced by substitution of (226) given that a r (x) P d (x) = 0 by definition to 

V x =-ye 2 (TL) (231) 

Note that (231) is guaranteed to be negative definite independent of the choice on 
the null-control vector y with proper selection of y. 

Given that the null-control vector y can be arbitrarily chosen, it can be selected to 
ensure V 2 <0 which, when combined with (231), ensures that V is negative definite and 
global asymptotic stability of the underactuated rigid spacecraft system is guaranteed. 
The corresponding time-derivative of V 2 is 


V 2 =k-'o)l6) a -2q 4 

= fc"V^ 2 i (®)°\ + k ^^l F n (®)®„ + + k~W a P d (x)y + q> 


(232) 
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Let the null-control vector y be defined as 


y = < ( x )[-<*>% (®H - «> - kq>] - r/oO fl (233) 

with a positive scalar gain d e M. Similar to the damping of the control coefficient 
generalized inverse, the damped actuated dynamics generalized inverse i\ d (x) given by 
(222) implies that 

lh;wi s 4- ( 234 > 

Pi 

and that 

lim L+(x)|| = 0 (235) 

(0 a — 

and that x\ + d (x) pointwise converges to r\ + d (x) as J3 2 vanishes. 

In the condition where to T a P d (x) = 0 lt2 , the null-control vector can be seen to not 

affect the stability of the system in a Lyapunov sense from (232). However, do to the 
arbitrarily small nature of this region by selection of the damping coefficient fP , it is 
conjectured that the system is globally stable to this arbitrarily small region defined by 
the damping coefficients /( and (P with the feedback control law 

u = a^(x)Z7(x) + P^(x)y (236) 

where y is given by (233). 
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E. FEEDBACK-LINEARIZING QUATERNION FEEDBACK REGULATOR 


Let the null-control vector y be chosen as 

y = -F 2I (to) co u - F 22 (to) o a - ck o fl - kq p (237) 

where d,kare positive scalar gains to be determined. In this manner, substitution of 
(237) into (212) yields the class of control laws 

u, =u d +P d (x)[-F 2l ((o)co u -F 22 ((o)(o a -d<o a -kq p ] (238) 

The close-loop actuated dynamical subsystem becomes 

<*>« = F 21 (to) oj u + F 22 (to) to„ + u, ; + P d (x)[-F 21 (to) - F 22 (ra)(o a - dto a - kq p ] (239) 

which converges by (210) and (214) to the linear dynamical system (Bajodah 2009a) 

™ a =- d ™ a - k % ( 24 °) 

Furthermore, the linear dynamical system given by (240) is upper bounded by 
selection of the damping coefficient /?, forming a domain of attraction in x whereby 

||a(x)||> ff due to the damped controls coefficient generalized inverse defined by (209). 

Therefore, the dynamical system given by (240) coupled with the remaining equations of 
motion given by (217) through (219) form a feedback linearizing transformation from the 
underactuated rigid spacecraft angular motion equations given by (175) over a domain of 
attraction in x. 

Conjecture 2. Let the function <fr(x u ) be globally twice continuously 
differentiable and satisfy the condition given by (184), the controls coefficient a(x) given 
by (190) and controls load b(x)given by (191) relative to ^(x M ) of the linear unactuated 
second-order dynamics given by (185) along the drift vector /'(x) = F(x)co ;i/ of the 
underactuated rigid spacecraft system of angular motion equations given by (175) 
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Furthermore, let the damped generalized inverse of the controls coefficient a), (x) he 
given by (209), the damped null-projector of the controls coefficient P d (x)be given by 
(212), and the null-control vector y be given by (221). If the zero-actuated state 
Jacobian J {) (x} given by (193) is nonsingular, then the positive scalar feedback gains 
k,d eR can be chosen to stabilize the closed-loop underactuated rigid spacecraft system 
of angular motion equations given by (217) through (220) by bounded, smooth, time- 
invariant state feedback control to an arbitrarily small vicinity of the equilibrium of the 
system defined by the damping coefficients /( within a domain of attraction in x(f). 

Proof: Assuming the inverse damping gain k 1 e R exists and is positive 
definite, consider the control Lyapunov function 


1 2 

1(0 a + 9i + 9 2 + c h + (<?4 - !) 


(241) 


Note that V is positive definite. This particular Lyapunov function is well known 
for spacecraft attitude dynamics appearing in Wie and Barba (1985); Wie, Weiss and 
Arapostathis (1989); and Cristi and Burl (1993) among others. 

Differentiating V along the trajectories of the underactuated system given by 
(175) yields 


V = (3) T a k x (i) a -2q A 
= -(O r a k 1 dGi a +q u (O u 


(241) can be rewritten as 


v=L|Kf+2(i-« 4 ) 


(242) 


(243) 


and (242) can be rewritten as 
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( 244 ) 


V = -kd\\(o a f +q u eo u 
Given ^(x M ) defined by (183) such that 

</>{\) = a >u+ a( lu ( 245 ) 

(244) can be rewritten as 

V = 7XA) -aql - kd ||co a f (246) 

From Lemma 1 in Appendix B, ||co a (t)|| is uniformly bounded, 
i.e,. |to a (t) | < c Vf. By defining /(/)el to be the scalar function described by 

f(t) = \\\(O a {r)fdT (247) 

0 

and given the designed dynamics in (185), ^(x i( )(r)^0 exponentially, application of 
Lemma 2 in Appendix B yields 

lim/(t) = /(oo)<+oo (248) 

and 

t 

/(t) = J|(n fl (r)||" Jr <c Vt (249) 

0 

Furthermore, with (247) and (249) defined, Barbalat’s Lemma (Slotine and Li 
1991) can be applied by conjecture given an arbitrarily small selection of the damping 
coefficient J3 1 on the controls coefficient generalized inverse to yield 

lim f (t)= lim ||(d (t)\\ = 0 (250) 

f ->+00 V 7 t —^+00 ' V 7 II 

From Lemma 3 in Appendix B, 
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( 251 ) 


SsIMOM 

and given Lemma 4 in Appendix B, (250) and (251) yield 

limllq (r)|| =° (252) 

oo II P V 7 II 

Thus, the closed-loop equations of motion given by (217) through (220) are stable 
within a domain of attraction in x(t). 

F. DESIGN EXAMPLE 

In order to demonstrate the effectiveness of the presented quaternion feedback 
regulators with both null-control vector selections for underactuated spacecraft, a control 
design example is presented that considers an asymmetric rigid spacecraft with principal 
moments of inertia (in kgm) of J n =32.5, J 22 = 25 and J 33 =12.5 . These selections of 

moments of inertia correspond to a 30 kg rigid rectangular spacecraft with uniform mass 
distribution and sides of length 1x2x3 m. For purposes of representing real-world 
spacecraft torqueing capabilities, the maximum available torque about each control axis 
is limited to 1 Nm. 

From (190), the control coefficient a(x) relative to the feedback linearizing 
transformation </>(x u ) given by (183) is 


l ( x ) = 


_i 

z J n 

Z J 11 


The determinant of the associated (x) as defined by (193) is 


(253) 
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which is non-zero for all x(7) e R 7a1 . 

For the numerical simulations that follow, a fourth-order Runge-Kutta numerical 
integration scheme is used to integrate (217) through (220) with a fixed time step of .Is. 
Appendix C provides the code for the algorithm implemented via an Embedded 
MATLAB® function. In order to accommodate exact zero initial conditions on to and 
q without failure of the regulator, the measured state is perturbed to 1x10 Mn these 
situations. 


1. Numerical Simulations using Lyapunov Function-Based Quaternion 
Feedback Regulator 

a. Rest to Rest Reorientation from Large Initial Angle Offset Given 
Low Available Control Torque 

A large-angle rest-to-rest reorientation maneuver is considered. The 
quaternion describing the initial orientation ( t 0 ) of 2f with respect to I is taken to be 

q(t 0 ) = [.57,.57,.57,.159] r which, by considering a 3-2-1 sequence of rotations, 

corresponds to initial Euler angles of [or, (t 0 ),a 2 (t 0 ),a 3 (t 0 )] = [109.8,-27.9,109.8]deg . 

The final orientation is taken to be q^ ) = [0,0,0, l] r . The initial and final angular 

velocity vector is co(t 0 ) = = [0,0,0] r . This reorientation maneuver is chosen in 

order to compare against the fully actuated quaternion feedback controller presented in 
Wie and Barba (1985). The chosen gains are 

y = .70, a = 1.25, d = 7.5, k = 2.25, f v f 2 = lxl (T 4 
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Figure 15 shows the time histories of the Euler angles and angular rates and 
demonstrates smooth attitude stabilization within 90 s. Figure 16 shows the time 
histories of the control torques. The control noise near the inherent singularity in the 
null-control vector given by (233) is evident throughout the maneuver until the 
stabilization region is reached. 
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Figure 15. Time Histories of the Euler Angles and Angular Velocities for a Rest to 
Rest Reorientation Manuever from Large Initial Angle Offset Using the Lyapunov- 
Function Based Control Law and Nominal Available Control Torque 
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Figure 16. Time Histories of the Control Torques for a Rest to Rest Reorientation 
Manuever from Large Initial Angle Offset Using the Lyapunov-Function Based Control 
Law and Nominal Available Control Torque 
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b. Large Initial Angular Rate to Rest Stabilization and 
Reorientation Given Low Available Control Torque 

A three-axis detumbling and reorientation maneuver is considered. The 
quaternions describing the initial and final orientation of S d with respect to I are 

q^) = q(? 0 ) = [0,0,0,l] r and the initial angular velocity vector is taken to be 

a(t 0 ) = [l,-l,l]rad/s 2 . The final angular velocity vector is a[t f ) = [0,0,0] r . By 

assuming the same choice of parameters given above, Figure 17 depicts the time histories 
of the Euler angles and angular rates and demonstrates attitude stabilization within 300 s. 
The failure to smoothly converge to the origin is due to the limited control torque. Figure 
18 shows the time histories of the control torques and again shows significant control 
effort noise during the maneuver until the stabilization region is reached. 
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Figure 17. Time Histories of the Euler Angles and Angular Velocities for a Large 
Initial Angular Rate to Rest Reorientation Manuever Using the Lyapunov-Function 
Based Control Law and Nominal Available Control Torque 
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Figure 18. Time Histories of the Control Torques for a Large Initial Angular Rate to 
Rest Reorientation Manuever Using the Lyapunov-Function Based Control Law and 

Nominal Available Control Torque 
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c. Large Initial Angular Rate to Rest Stabilization and 
Reorientation Given Large Available Control Torque 

A three-axis detumbling and reorientation maneuver is considered again 
with increased torque capabilities as is typically studied in most literature involving 
underactuated spacecraft control. The quaternions describing the initial and final 

orientation of 2f d with respect to 7 are again taken to be = q(r 0 ) = [0,0,0,l] r and 

the initial angular velocity vector is again taken to be co(r 0 ) = [1,-1,l] rad / s 2 . The final 
angular velocity vector is to (t/) = [0,0,0] r . The maximum torque is taken to equal to the 

moment of inertia about the second principal axis so that u 2 is limited to 1 rad Is 2 . By 

assuming the same choice of parameters given above, Figure 19 shows the time histories 
of the Euler angles and angular rates, demonstrating attitude stabilization in less than 90s. 
Figure 20 shows the time histories of the control torques. The control noise is clearly 
evident near the inherent singularities in the null-control vector and its chattering effect 
on the angular rates. 
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Figure 19. Time Histories of the Euler Angles and Angular Velocities for a Large 
Initial Angular Rate to Rest Reorientation Manuever Using the Lyapunov-Function 
Based Control Law and Large Available Control Torque 
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Figure 20. Time Histories of the Control Torques for a Large Initial Angular Rate to 
Rest Reorientation Manuever Using the Lyapunov-Function Based Control Law and 

Large Available Control Torque 


109 

















































d. Rest to Rest Reorientation Manuever Given an Initial Angle 
Offset about Only the Unactuated Axis and Low Available 
Control Torque 

An attitude error stabilization maneuver is considered about the 
unactuated body-fixed axis with low torque capabilities. The quaternion describing the 

initial orientation (t 0 )of with respect to I is taken to be q(f 0 ) = [.17,0,0,.98] T which, 
by considering a 3-2-1 sequence of rotations, corresponds to initial Euler angles of 
[eq (t 0 ),oc 2 (t 0 )’ a 3 (t 0 )] = [20,0,0]deg. The final orientation is taken to be 

q [t f j = [0.0.0,lf. The initial and final angular velocity vector is 

co(t 0 ) = = [0,0,0] r . By assuming the same choice of parameters given above, 

Figure 21 shows the time histories of the Euler angles and angular rates, demonstrating 
attitude stabilization in less than 70 s. Figure 22 shows the time histories of the control 
torques and demonstrates significant control noise near the singularities in the control 
algorithm despite the small initial angular offset. 
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Figure 21. Time Histories of the Euler Angles and Angular Velocities for a Rest to 
Rest Reorientation Manuever Given an Initial Angle Offset about Only the Unactuated 
Axis Using the Lyapunov-Function Based Control Law and Nominal Available Control 

Torque 
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Figure 22. Time Histories of Control Torques for a Rest-To-Rest Reorientation 
Manuever Given an Initial Angle Offset about Only the Unactuated Axis Using the 
Lyapunov-Function Based Control Law and Nominal Available Control Torque 
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e. Attitude Maintenance in the Presence of Relatively Large 
Disturbance Torques 

In order to demonstrate a limitation of the presented Lyapunov-function 
based quaternion feedback regulator with respect to attitude maintenance, an attitude 
error maintenance scenario is considered with low torque capabilities in the presence of 
relatively large disturbance torques. The quaternion describing the initial orientation 

(t 0 ) of 3B d with respect to I is taken to be q (f 0 ) = [0,0,0, l] ; while the final orientation is 
taken to be ) = [0,0,0, l] r . The initial and final angular velocity vectors are 

co(t 0 ) = = [0,0,0] r . The same choices of parameters given above are used to 

include a maximum torque about the two control axes of 1 Nm. Furthermore, a relatively 
large disturbance torque of .1 Nm is imparted about each body axis. For the presented 
rigid spacecraft inertia matrix of / = diag (32.5,25,12.5) kgm 2 , this corresponds to 

angular accelerations of .003m/s 2 , .004m/s 2 and .008 m/s 2 about the respective body 
axes. Figure 23 shows the Euler angles and angular rates and demonstrates that although 
the desired attitude is unable to be maintained, it does remain bounded. Figure 24 reports 
the control history for the attitude maintenance scenario and again shows the large 
control effort noise given the near-singularities that exist in the controller. 
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Figure 23. Time Flistories of Euler Angles and Angular Velocities for an Attitude 
Maintenance Scenario in the Presence of Large Disturbance Torques Using the 
Lyapunov-Function Based Control Law and Nominal Available Control Torque 
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Figure 24. Time Histories of the Control Torques for an Attitude Maintenance 
Scenario in the Presence of Large Disturbance Torques Using the Lyapunov-Function 
Based Control Law and Nominal Available Control Torque 
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2. Numerical Simulations using Feedback Linearized Based Quaternion 

Feedback Regulator 

a. Rest to Rest Reorientation from Large Initial Angle Offset and 
Low Available Control Torque 

A large-angle rest-to-rest reorientation maneuver is considered. The initial 
and final conditions of the spacecraft are taken to be equivalent to the considered 
maneuver above with the Lyapunov-based approach. The chosen gains are 

Y = -70, a = 1.25, d = 7.5, k = 2.25, /?, = lxl O' 4 

Figure 25 shows the time histories of the Euler angles and angular rates, 
demonstrating attitude stabilization within 100 s. The discontinuity in the Euler angles is 
due to singularities in the Euler rotation sequence and does not affect the algorithm. 
Figure 26 shows the time histories of the control torques and demonstrates a significant 
advantage over the Lyapunov-function based controller with no control noise. 
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Figure 25. Time Histories of the Euler Angles and Angular Velocities for a Rest to 
Rest Reorientation Manuever from Large Initial Angle Offset Using the Feedback 
Linearizing Control Law and Nominal Available Control Torque 
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Figure 26. Time Histories of the Euler Angles and Angular Velocities for a Rest to 
Rest Reorientation Manuever from Large Initial Angle Offset Using the Feedback 
Linearizing Control Law and Nominal Available Control Torque 
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b. 


Rate Stabilization and Reorientation (Low Torque) 


A three-axis detumbling and reorientation maneuver is considered. The 
initial and final conditions of the spacecraft are taken to be equivalent to the considered 
maneuver above with the Lyapunov-based approach. By assuming the same choice of 
parameters, Figure 27 depicts the time histories of the Euler angles and angular rates, 
demonstrating convergence to the zero state in less than 250 s. Figure 28 shows the time 
histories of the control torques and again shows significant control effort savings over the 
Lyapunov-based approach without sacrificing realizability for the system. 
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Figure 27. Time Histories of the Euler Angles and Angular Velocities for a Large 
Initial Angular Rate to Rest Reorientation Manuever Using the Feedback Linearing 
Control Law and Nominal Available Control Torque 
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Figure 28. Time Histories of the Control Torques for a Large Initial Angular Rate to 
Rest Reorientation Manuever Using the Feedback Linearizing Control Law and Nominal 

Available Control Torque 
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c. 


Rate Stabilization and Reorientation (Large Torque) 


A three-axis detumbling and reorientation maneuver is considered again 
with increased torque capabilities as is discussed for the Lyapunov-based approach. The 
initial and final conditions of the spacecraft are taken to be equivalent to the considered 
maneuver above with the Lyapunov-based approach. By assuming the same choice of 
parameters given above to include a maximum torque about the second and third body 
axes so that u 2 is limited to 1 rad / s 2 , Figure 29 shows the time histories of the Euler 

angles and angular rates. Due to the large available torque, the spacecraft is stabilized to 
the desired attitude in less than 100 s. Figure 30 shows the time histories of the control 
torques and shows that even with a large available torque, the control does not saturate 
during the maneuver as is the case with the Lyapunov-based solution. 
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Figure 29. Time Histories of the Euler Angles and Angular Velocities for a Large 
Initial Angular Rate to Rest Reorientation Manuever Using the Feedback Linearizing 
Control Law and Large Available Control Torque 
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Figure 30. Time Histories of the Control Torques for a Large Initial Angular Rate to 
Rest Reorientation Manuever Using the Feedback Linearizing Control Law and Large 

Available Control Torque 
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d. Unactuated Attitude Error Stabilization (Low Torque) 

An attitude error stabilization maneuver is considered about the 


unactuated body-fixed axis with low torque capabilities. The initial and final conditions 
of the spacecraft are taken to be equivalent to the considered maneuver above with the 
Lyapunov-based approach. By assuming the same choice of parameters given above, 
Figure 31 shows the time histories of the Euler angles and angular rates, demonstrating 
smooth reorientation to the desired attitude in less than 80 s. Figure 32 shows the time 
histories of the control torques and shows the advantage in terms of control effort over 
the Lyapunov-based solution. 
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Figure 31. Time Histories of the Euler Angles and Angular Velocities for a Rest to 
Rest Reorientation Manuever Given an Initial Angle Offset About Only the Unactuated 
Axis Using the Feedback Linearizing Control Law and Nominal Available Control 

Torque 
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Figure 32. Time Histories of the Control Torques for a Rest to Rest Reorientation 
Manuever Given an Initial Angle Offset About Only the Unactuated Axis Using the 
Feedback Linearizing Control Law and Nominal Available Control Torque 
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e. Attitude Maintenance in the Presence of Relatively Large 
Disturbance Torques 

An attitude error maintenance scenario is considered with low torque 
capabilities in the presence of relatively large disturbance torques using the presented 
feedback linearizing control law. By assuming the same initial and final conditions and 
parameters used for the Lyapunov-based approach. Figure 33 shows the Euler angles and 
angular rates and demonstrates that although the desired attitude is unable to be 
maintained, it remains bounded. By comparison to Figure 23 for the Lyapunov-function 
based control law, the feedback linearizing control law provides an additional measure of 
stability by affecting a strictly periodic nature to the spacecraft’s attitude and a tighter 
bounds on the attitude error departures. Figure 34 reports the control history for the 
attitude maintenance scenario and depicts a relatively smooth control effort throughout 
the maneuver. 
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Figure 33. Time Histories of Euler Angles and Angular Velocities for an Attitude 
Maintenance Scenario in the Presence of Large Disturbance Torques Using the Feedback 
Linearizing Control Law and Nominal Available Control Torque 
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Figure 34. Time Histories of the Control Torques for an Attitude Maintenance 
Scenario in the Presence of Large Disturbance Torques Using the Feedback Linearizing 
Control Law and Nominal Available Control Torque 
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y. 


CONCLUSION 


A. SUMMARY OF CONTRIBUTIONS 

A planar laboratory testbed is first introduced for the simulation of autonomous 
proximity maneuvers of uniquely control actuator configured spacecraft. This testbed 
consists of a floating robotic simulator via planar air bearings on a flat floor that is 
equipped with dual vectorable cold-gas thrusters and miniature control moment 
gyroscope. Inertial position and attitude measurements are obtained and fused with a 
discrete Kalman filter and linear quadratic estimator for navigation while feedback 
linearized control coupled with a linear quadratic regulator and Schmitt trigger logic 
directly command the control moment gyro and command dual in-plane vectorable cold- 
gas thrusters through a Pulse Width Modulation scheme. 

The presented experimental tests of autonomous closed path proximity maneuvers 
of the spacecraft simulator offer significant sample cases. The experimental results, 
which show good repeatability and robustness against disturbance and sensor noise, 
validate the proposed estimation and control approaches and validate, in particular, the 
analytical small-time local controllability of the uniquely configured system. The 
achieved accuracy in following the reference trajectory (respectively, ~ 1 cm for 
translation and ~ .5 deg for rotation given only the vectorable thrusters as control inputs) 
demonstrates both a feasible and promising actuator configuration for small spacecraft 
design given size, weight and fuel storage considerations. 

Leveraging the research on the 3-DoF robotic testbed, a rigorous small-time local 
controllability analysis is conducted to determine the minimum number of control 
actuators for full-order spacecraft in a proximity operations environment. Two potential 
configurations flow from this analysis. The first involves a traditional approach of fixed 
thrusters whereby two fixed thrusters are oppositely mounted and used to translate the 
spacecraft while two sets of paired thrusters provide attitude control. The second 
configuration involves a novel design consisting of only two hemispherically vectorable 
thrusters. By combining the vectorable thrusters in an opposing manner, the typical 
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spacecraft system requiring six or more control actuators to achieve 6-DoF control is 
reduced to only two. Furthermore, by simply adding two additional paired thrusters to 
provide a torque about the vector that joins the two vectorable thrusters, this system 
becomes feedback linearizable and similar methods to those presented for the 3-DoF 
robotic testbed can be employed to control both its attitude and position. 

By specifically focusing on the minimal number of control actuators for the 6- 
DoF system, that is to say where only two vectorable thrusters are required to provide 
small-time local controllability, a control problem exists where one must achieve three- 
axis stabilization of the attitude with only two control torques. To solve this problem, a 
novel quaternion feedback regulator is presented that capitalizes on recent developments 
in generalized inversion and perturbed feedback linearizing control. A desired second- 
order linear dynamics in a function of the angular velocity and error quaternion 
components about the unactuated axis is evaluated along trajectories of the 3-DoF 
spacecraft angular motion equations that include both Euler’s dynamic equations and the 
error quaternion kinematics. The control variables, which are composed of a particular 
and auxiliary part, are determined to yield a stable underactuated system. The particular 
part seeks to realize the desired linear dynamics through generalized inversion of the 
controls coefficient while the null-control vector in the auxiliary part is selected to yield 
stability of the full underactuated system. Two control design methodologies are derived 
that involve separate constructions of the null-control vector. The first is Lyapunov 
based and yields a stable underactuated spacecraft system while the second involves 
perturbed feedback linearization and yields stability of the system within a domain of 
attraction. For both cases, in order to overcome the potential for the control laws to 
produce numerical instability as the generalized inverse of the controls coefficient 
becomes singular, the generalized inverse is damped near a bound on the singularity, thus 
providing smooth control and resulting in stabilization near to the zero error state. As 
such, neither controller contradicts the conjecture by Byrnes and Isidori (1991) that no 
smooth, time-invariant state feedback controller can be found to locally asymptotically 
drive the system to the zero error state. However, given the arbitrary selection of 
damping coefficients, the region about the zero error state can be reduced to an 
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insignificant sized ball. Several simulations results are presented for an underactuated 
asymmetric spacecraft with two, realistically bounded, body-fixed torques to demonstrate 
the wide-reaching capabilities of both presented null-control vector designs. Included are 
a large-angle rest-to-rest reorientation maneuver, a large angular velocity detumbling and 
simultaneous reorientation maneuver and a single-axis reorientation about the unactuated 
axis. The proposed attitude control method is not intended to provide attitude 
maintenance or for attitude tracking in the presence of relatively large disturbance 
torques; however, it may prove widely applicable to detumbling and reorientation 
maneuvers of spacecraft with only two available control torques. 

B. POSSIBLE FUTURE DEVELOPMENTS 

The following iterating research remains open: 

1. Validate by numerical and/or experimental methods the feedback linearizing 
control method proposed in Chapter III for the full 6-DoF relative dynamics 
problem. 

2. Determine a roto-translation control strategy for a spacecraft with two 
oppositely mounted hemispherical thrusters. A possible solution is to 
investigate using a Frenet frame waypoint tracking solution where the forward 
face is first made normal to the path by one of the quaternion feedback 
regulators presented in Chapter IV and then a separate controller is used to 
translate the spacecraft. 

3. Provide a formal proof of global stability for both the Lyapunov-function 
based and feedback linearizing underactuated quaternion feedback controllers 
presented in Chapter IV. Given the presented numerical simulation results 
that demonstrate a very large domain of attraction, it is possible that global 
stability for both controllers can be rigorously demonstrated. 

4. Compare both underactuated attitude controllers to the optimal solution for 
each presented case. In this manner, a method of gain tuning may be achieved 
to ensure a near-optimal solution can be achieved. 
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APPENDIX A. MATLAB CODE TO PERFORM SMALL TIME 
LOCAL CONTROLLABILITY STUDIES 


function [L, CO] = nlctrb(fx,Gx,P) 


%NLCTRB Compute the controllability matrix and associated 

% Lie algebra for a nonlinear control-affine system 

% (with or without drift) of the following form: 

'o 

% xdot = f(x) + G(x)u, f(P) = 0, x->RNx, u->RNu 

"6 

% The following three cases pertain for the condition of the 
% dimensions of the Lie Algebra Matrix L where Nu is 

% the dimension of the control space and Nx is the 

% dimension of the state space: 


1. dim(L(Delta)) = Nu 

2. Nu < dim(L(Delta)) < Nx 

3. dim(L(Delta)) = Nx 


The system is completely integrable 
The system is nonholonomic, but not 
Small-Time Locally Controllable (STLC) 
The system is nonholonomic and STLC 


if the system is driftless (i.e. f_0(x) 
the zero vector for fO(x). 


= 0 for all x in X) enter 


Evaluate this script to establish symbolic objects for state 

variables prior to calling function where state = number of 
state variables 

syms X 

for a = 1: state 

x = strcat('x',num2str(a)); 
syms(x); 

X (a) = x; 

end 


This algorithm uses the P.Hall Basis algorithm as explained in 
LaValle, S.M. Planning Algorithms, 2006. 


%Developed by Jason S. Hall 
%Naval Postgraduate School 
%Feb 2007 

%Modified July 2009 to work with new Symbolic Matlab Toolbox 

oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

%% State Variable Assignment 
[Nx,Nu] = size(Gx); 

X = sym(zeros(Nx, 1)) ; 

for a = 1:Nx 

x = strcat( 1 x ’,num2str (a)); 
syms(x); 

X (a) = x; 

end 

%% Vector Variable Assignment 
syms h fO 
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1 


if isequal(fx,zeros(Nx,1)) == 
dimh = Nu; 


q {l} 

= 'hi'; 


V 

= genvarname('h' 

,who); 

eval([v 

' = Gx(:,l) ’]) 


L {1} 

= q{ i} ; 


for a = 

2:dimh 


q{ a ] 

\ = strcat('h' 

,num2str(a)); 

V 

= genvarname('h',who); 

eval([v ' = Gx(:,a); 

' ]); 

L(a) 

= q (a) ; 


end 



st 

= i; 


dimh 

= Nu + 1; 


q {1} 

= 'fO ' ; 


fO 

= fx; 


st 

= 0; 


for a = 

2:dimh 



q{a} = strcat( 'h' ,num2str(a-1)) 
v = genvarname( 'h' ,who) ; 

eval([v ' = Gx(:,a-1); '] ) ; 

L(a-1) = q(a) ; 

end 

end 

dim = ones(1,Nu+1); 

Xn = zeros(Nx,1); 

for a = 1:Nx 


if P(a) ~= 0; 



Xn(a) = randn; 

end 


end 


CO 

= Gx; 

Ldel_x 

= double(subs(Gx,X,Xn)); 

rhmat 

= rank(Ldel_x); 

%% Required Index Assignment 

i—1 
i—1 

M 

r* 

= size(Gx); 

k 

= 1; 

j j 

= i; 

%% Algorithm 

for a = 

2:Nu+5 


[o, p] = size (q) ; 
if rhmat == Nx 


break 

end 

for b = st:Nu 

if rhmat == Nx 
break 

end 

for c = 1:p 

if rhmat == Nx 
break 

end 

if b == 0 
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rl = 'f0' ; 

else 

rl = strcat( 'h' ,num2str(b)); 

end 

if iscellstr(q(a-1,c)) && isequal(char(q(a-1,c)),' n 
r2 = strcat(q(a-1,c)); 
r = strcat(rl,r2); 
s = sym(char(r)); 

else 

r = { ' n ' }; 
s = sym(char (r)); 

end 

if strcmp(rl,r2) 
r = { ' n ' } ; 
s = sym(char(r)); 

elseif (strcmp(rl, 'hi' ) && strcmp(r2, 'f0, hi' )) I ! . . 

(strcmp(rl, 'h2' ) && strcmp(r2, 'f0,h2' )) 

r = { ' n ' }; 
s = sym(char (r)); 

else 

if size(rl,2) <= 3 
r3 = char(r2); 

if str2double(rl(2:end)) >= ... 

str2double(r3(2:end)) && ... 
isequal(r3(end) ,' n' ) == 0 
r = { ' n ' }; 
s = sym(char(r)); 

end 

end 

if isequal(r, {'n '})~=1 
q (a, k) = r; 
k = k + 1; 

end 

if length(rl) >= 2 && isequal(char(s) ,' n' ) ~= 1 

s2 = s; 

end 

for d = 1:dimh 

if size(rl) == size(char(q(1,d))) 
if rl == char(q(l,d)) 
for e = 1:Nu+1 

if length(s) == e+1 

si {e } (dim (e) , : ) = s; 
kk = dim(e); 
dim(e) = dim(e) + 1; 

end 

end 

end 

end 

end 

end 

if exist( 's2' , 'var' ) 
fng = eval(s2); 
zz = size(fng,2); 
while zz > 1 

zz = zz -1; 
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fv = fng(:,zz); 
g = fng(:,zz+1); 
dfdx = jacobian(fv,X); 
dgdx = jacobian(g,X); 
fng(:,zz) = dgdx*fv-dfdx*g; 

end 

brack = fng(:,zz); 

FNG = vpa(subs(brack,X,Xn),2); 

FNG1 = subs (brack, ' X' , 'P') ; 

hmat = [Ldel_x,FNG]; 

pp = size(hmat,2); 

rhmat = double(rank(hmat)); 

if rhmat == pp 

Ldel_x = hmat; 


[CO,FNG1]; 

char(si{end}(kk,:)); 
cellstr(qn(9:size(qn,2)-2)) 
11 + 1 ; 


CO 

qn 


L(11+1) 
11 


end 

if rhmat == Nx 
break 

end 

clear s2 


end 

end 

end 

jj = jj+i; 

k = 1; 

end 

PHall_Basis = o; 

if Nu < Nx && isequal(fx,zeros(Nx,1)) == 1 
if rhmat == Nx 

fprintf('The system is STLC \n \n') 
elseif rhmat == Nu 

fprintf('The system is completely integrable \n \n') 
elseif rhmat < Nx && rhmat > Nu 

fprintf('The system is not STLC \n \n') 

else 

fprintf('The system doesn''t fit \n \n') 

end 

elseif Nu < Nx && isequal(fx,zeros(Nx,1)) == 0 
if rhmat == Nx 

fprintf('The system is at least accessible, ') 
fprintf (' verify no bad brackets to determine STLC\n \n') 
elseif rhmat == Nu 

fprintf ('The system is completely integrable \n \n') 
elseif rhmat < Nx && rhmat > Nu 

fprintf ('The system is not accessible \n \n') 

else 

fprintf ('The system doesn''t fit \n \n') 

end 

end 

fprintf(' L(Delta) = span {%s' ,char(L(1))) 
for a = 2:length(L) 
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fprintf(' ,%s' ,char(L(a))) 

end 

fprintf('} \n' ) 

fprintf(' \nThe LRAC is %i \n', rhmat); 

fprintf(' P.Hall Basis to Depth %i \n' ,PHall_Basis) 
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APPENDIX B. PROOFS OF EMPLOYED LEMMAS 


Lemma 1 . (o a (t) is uniformly bounded, i.e. to a (t) <c\/t 


1 m2 

Proof. By rearranging V = — k |( 0 a |[ + 2(l - q 4 ) to yield 


£||(d u || 2 =2V-4(1-4 4 ) 


(255) 


and then substituting this into V = q h (j ){\ u ) - aq] - kd ||a> a || 2 produces 

V = ~ 2dv + 4d(l-q 4 ) + qj(x u ) - aq 2 u (256) 

The desired stable linear dynamics ^(x H ) + 2;y0(x M ) + 7 2 0(x i( ) = 0 and the 
condition det[y o (x)]^0 for the infinite set of control laws u rf = u rf +P d (x)y with 
u rf =a + d (x)&(x) ensure that 

lim^fx ) = 0 (257) 

t —^oo v 7 


for an arbitrarily small selection of that damping coefficient f\ on the controls coefficient 
generalized inverse a^(x). Therefore, given c/ 4 | < 1 Vf from the definition of the error 
quaternion (Cristi and Burl 1993), V (7) is bounded and |e> a (t)|| is bounded. 

Lemma 2. If 0(x m )(7)—>0 exponentially, i.e. |^(x M )(t)| < Aet > 0 for some A, y > 0, 
then 

lim/(0 = /(oo)<+oo (258) 

f->+oo V 7 V 7 

where f (t) is defined as 
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(259) 


/M=jlkMlf dr 

o 

Proof. For all t e R (Cristi and Burl 1993) 

t 

y(r) = y( 0) + jV(r)dr (260) 

0 


S 2 

Substitution of the Lyapunov function V given by V = — fc|<o a || + 2(l-g 4 ) and its 
derivative V = q u </)(* u ) - aq 2 u - kd ||co a || 2 yields 

y (0 = y (°)“ M lll co «( r )t d* + \q le { T )<P{ x u){ T )dT-a\q;, (r)dr (261) 

0 0 0 


Since V (?) > 0 V? 


kd‘\\\(o a (r)||" dz < V (0) + j|^(x„ )(r )|dr (262) 

0 0 


Letting 2 = 1 /kd yields 


Ilk Mil 2dz - AV (°) +A l k( x »)Mk 


o 


o 


(263) 


and since ^(x H ) exponentially decays through the selection of y in 
^(x i( ) + 2y^(x i( ) + 7 2 ^(x u ) = 0 to satisfy the characteristic equation given 
by s 2 + 2 ys + y 2 = 0, the integral is finite. Then for some constant c 

t 

/(V) = J||(D fl (r)|~ Jr <c Vf (264) 

0 

Given (264) is monotonically non-decreasing and upper bounded, it has a finite 
limit as t —> +oo and thus for some constant c 
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-I-UU 

/(°o)= J|to fl (r)| 2 dr<c 


(265) 


Lemma 3. Given (257), 


lim IK x «)ll =G >u +a <i« =° 

t —^co v 7 


(266) 


which implies that 


lim k ( 0 | = -^„(0 


(267) 


Furthermore, with the lim to (t) =0 from (250) and the kinematic equation for 

t—> oo V / ■ 


the scalar component of the error quaternion given by q A - 0 x u + Hp^a 


® 2 


(268) 


is implied. Therefore, substitution of (267) into (268) yields 


lim|4 4 | = --C 

f->co 2 


(269) 


so that 


lim qAt ) =0 

/—»co ' 7 


(270) 


Additionally, given the kinematic equation for the unactuated component of the 
error quaternion expressed by q u = fq 4 co u - q j ,the lim |a> a (/ )|| = 0 from (250) and 
the result from (270) yields 


lim du (0 = o 


(271) 


And, therefore, from (267), 
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(272) 


lim <» (t)\ = 0 

t —^oo V 7 1 

and furthermore 

limllx (t)\\ = 0 (273) 

?-»oo '• V 7 H 

Lemma 4. Given limllco (t)\\ = 0 from (250) and limllx (7)11 = 0 from (273), and given 
the kinematic equations for the error quaternion as expressed by 

*!/> = ®2xl ^21 ) X h + ^2lifi.B d B )®a 

limllq (r)|| =0 (274) 

Proof. From q p =[0 2jtl F 2l (q B<fB )]x u + F 22 (q^ B )o> fl , given lim|o a (f)|| = 0 

from (250) and limllx (f)|| = Ofrom (273) 

t —^oo V 7 H 

Hk(f)|=0 (275) 

f—>oo 11 1 v II 

(b a = -dt o fl -kq p can be rewritten as 

(™ a +V<lp) + l <( G >a+V<lr) = V<lp ( 276 ) 

where k = d > 0, o = kd l > 0. By defining 

e M=®«M +a MO ( 277) 

and substituting this into (276) yields 

e(r) = -/re(r) + crq ;; (r) (278) 

From (275) and given k > 0 

lim|e(0|l = 0 (279) 

oo H V 7 H 
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and, therefore, substitution of (279) into (277) yields 


lim q (t) =0 


(280) 
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APPENDIX C. MATLAB CODE FOR UNDERACTUATED 
QUATERNION FEEDBACK CONTROLLER 


function [T] = RateStable(w, q, qd, wd, J, d, k, gamma,... 

Betal, Beta2, al, select) 

% Quaternion Feedback Attitude Regulator for Underactuated Spacecraft 


for 


end 

for 


i = 2:4 
if q(i) == 0 

q(i) = Betal; 

end 

i = 1:3 

if w(i) == 0; 

w(i) = Betal; 

end 


end 





qO 

= 

q(i); 



qi 

= 

q (2) ; 



q2 

= 

q (3) ; 



q3 

= 

q (4) ; 



qdO 

= 

qd(1); 



qdl 

= 

qd(2); 



qd2 

= 

qd(3); 



qd3 

= 

qd(4); 



qdv 

= 

[qdl;qd2;qd3]; 



qv 

= 

[ql; q2; q3] ; 



qx 

= 

[0 -q3 q2;q3 0 -ql;-q2 ql 

0] ; 


qOe 

= 

q0*qd0 + qv'*qdv; 



qve 

= 

qd0*qv - q0*qdv + qx*qdv; 



qle 

= 

qve (1); 



q2e 

= 

qve(2); 



q3e 

= 

qve (3); 



we 

= 

w - wd; 



wle 

= 

we (1) ; 



w2e 

= 

we (2); 



w3e 

= 

we (3); 



wx 

= 

[0 -w3e w2e;w3e 0 -wle;-w2e wle 

0 

Jll 

= 

J(l,1);J12 = J(1,2);J13 = 

J (1,3 

); 

J21 

= 

J(2, 1);J22 = J(2,2);J23 = 

J (2,3 

); 

J31 

= 

J(3, 1);J32 = J(3, 2);J33 = 

J (3, 3 

); 

a 

= 

[1/2* (-2*wle*Jll*J13*J22+2 

*wle* 

Jl 


2*wle*J31*J22*J33-2*J12*w3e*J13*J32-2*J22*w3e*J23*J32-. 
2*J22*wle*J12*J23-2*J32*wle*J12*J33+2*w3e*J13*J12*J23+. 
2*w3e*J33*J23*J32+2*w3e*J12 A 2*J33+2*w3e*J22 A 2*J33+... 
2*wle*J13*J22 A 2+2*wle*J13*J32 A 2-2*w3e*J13 A 2*J22-. .. 
2*w3e*J22*J33 A 2+4*w2e*J12 A 2*J23+4*w2e*J23*J32 A 2-. .. 
4*w2e*J12*J13*J22-4*w2e*J32*J22*J33+. . . 
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2*wle*J31*J23*J32-al*q3e*Jll*J22*J33+al*q3e*Jll*J23*J32+. 
al*q3e*J21*J12*J33-al*q3e*J21*J13*J32-al*q3e*J31*J12*J23+ 
al*q3e*J31*J13*J22)/(Jll*J22*J33-J11*J23*J32-J21*J12*J33+ 
J21*J13*J32+J31*J12*J23-J31*J13*J22) ; 

1/2* (2*wle*Jll*J12*J33-2*wle*Jll*J13*J32+ . . . 
2*wle*J21*J22*J33-2*wle*J21*J23*J32+2*w2e*J12 A 2*J33+... 
2*w2e*J22 A 2*J33-2*w2e*J13 A 2*J22-2*wle*J12*J23 A 2-. .. 
2*w2e*J22*J33 A 2-4*w3e*J13 A 2*J32-4*w3e*J23 A 2*J32-. .. 
2*w2e*J12*J13*J32-2*w2e*J22*J2 3*J32+2*J13*w2e*J12*J2 3+. . . 
2*J23*wle*J13*J22+2*J33*w2e*J23*J32-2*wle*J12*J33 A 2+... 
2*J33*wle*J13*J32+4*w3e*J13*J12*J33+4*w3e*J23*J22*J33+... 
al*q2e*Jll*J22*J33-al*q2e*Jll*J23*J32-al*q2e*J21*J12*J33+ 
al*q2e*J21*J13*J32+al*q2e*J31*J12*J23-. .. 
al*q2e*J31*J13*J22)/(Jll*J22*J33-J11*J23*J32- . . . 
J21*J12*J33+J21*J13*J32+J31*J12*J23-J31*J13*J22) ] ; 

Lf2phi = l/2*al*wle*(-l/2*qle*wle-l/2*q2e*w2e-l/2*q3e*w3e)+... 

l/2*al*w3e*(l/2*q3e*wle+l/2*q0e*w2e-l/2*qle*w3e)- ... 
l/2*al*w2e*(-l/2*q2e*wle+l/2*qle*w2e+l/2*q0e*w3e)+... 

(1/2*al*q0e+(-(J12*J2 3-J13*J22)/(Jll*J22*J33- . . . 
J11*J23*J32-J21*J12*J33+J21*J13*J32+J31*J12*J23-. .. 
J31*J13*J22)*J21-(J12*J33-J13*J32)/(Jll*J22*J33- ... 
J11*J23*J32-J21*J12*J33+J21*J13*J32+J31*J12*J23-. .. 

J31*J13*J22)*J31)*wle+((J12*J33-J13*J32)/(Jll*J22*J33- .. 
J11*J23*J32-J21*J12*J33+J21*J13*J32+J31*J12*J23-. .. 

J31*J13*J22)*w3e+(J12*J23-J13*J22)/(Jll*J22*J33- . . . 
J11*J23*J32-J21*J12*J33+J21*J13*J32+J31*J12*J23-. .. 

J31*J13*J22)*w2e)*J11+((J22*J33-J23*J32)/(Jll*J22*J33- .. 
J11*J2 3*J32-J21*J12*J33+J21*J13*J32+J31*J12*J2 3-. . . 

J31*J13*J22)*w3e-(J12*J23-J13*J22)/(Jll*J22*J33- . . . 
J11*J23*J32-J21*J12*J33+J21*J13*J32+J31*J12*J23-. .. 

J31*J13*J22)*wle)*J21+(-(J22*J33-J23*J32)/(Jll*J22*J33- . 
J11*J23*J32-J21*J12*J33+J21*J13*J32+J31*J12*J23-. .. 

J31*J13*J22)*w2e-(J12*J33-J13*J32)/(Jll*J22*J33- ... 
J11*J23*J32-J21*J12*J33+J21*J13*J32+J31*J12*J23-. .. 

J31*J13*J22)*wle)*J31+(-(J12*J23-J13*J22 ) /(Jll*J22*J33- . 
J11*J23*J32-J21*J12*J33+J21*J13*J32+J31*J12*J23-. .. 

J31*J13*J22)*J22-(J12*J33-J13*J32)/(Jll*J22*J33- ... 
J11*J23*J32-J21*J12*J33+J21*J13*J32+J31*J12*J23-. .. 

J31*J13*J22)*J32)*w2e+(-(J12*J23-J13*J22)/(Jll*J22*J33- . 
J11*J23*J32-J21*J12*J33+J21*J13*J32+. .. 
J31*J12*J23-J31*J13*J22)*J23-(J12*J33-J13*J32)/ ... 

(J11*J22*J33-J11*J23*J32-J21*J12*J33+J21*J13*J32+... 
J31*J12*J23-J31*J13*J22)*J33)*w3e)*((((J12*J33-J13*J32)/ 
(J11*J22*J33-J11*J23*J32-J21*J12*J33+J21*J13*J32+... 
J31*J12*J23-J31*J13*J22)*w3e+(J12*J23-J13*J22)/ ... 

(J11*J22*J33-J11*J23*J32-J21*J12*J33+J21*J13*J32+. . . 
J31*J12*J23-J31*J13*J22)*w2e)*J11+((J22*J33-J23*J32)/ ... 
(J11*J22*J33-J11*J23*J32-J21*J12*J33+J21*J13*J32+... 
J31*J12*J23-J31*J13*J22)*w3e-(J12*J23-J13*J22)/ ... 

(J11*J22*J33-J11*J23*J32-J21*J12*J33+J21*J13*J32+... 
J31*J12*J23-J31*J13*J22)*wle)*J21+(-(J22*J33-J23*J32) / . . 
(J11*J22*J33-J11*J23*J32-J21*J12*J33+J21*J13*J32+. . . 
J31*J12*J23-J31*J13*J22)*w2e-(J12*J33-J13*J32)/ ... 

(J11*J22*J33-J11*J23*J32-J21*J12*J33+J21*J13*J32+. . . 
J31*J12*J23-J31*J13*J22)*wle)*J31) *wle+. .. 
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(((J12*J33-J13*J32)/(J11*J22*J33-J11*J23*J32-. .. 
J21*J12*J33+J21*J13*J32+J31*J12*J23-J31*J13*J22)*w3e+ ... 
(J12*J23-J13*J22)/(J11*J22*J33-J11*J23*J32-J21*J12*J33+ . 
J21*J13*J32+J31*J12*J23-J31*J13*J22)*w2e) *J12+. . . 

((J22*J33-J23*J32)/(Jll*J22*J33-J11*J23*J32-J21*J12*J33 + 
J21*J13*J32+J31*J12*J23-J31*J13*J22) *w3e-. .. 

(J12*J23-J13*J22)/(J11*J22*J33-J11*J23*J32-J21*J12*J33+ . 
J21*J13*J32+J31*J12*J23-J31*J13*J22)*wle) *J22+. .. 

(-(J22*J33-J23*J32)/(Jll*J22*J33-J11*J23*J32- ... 
J21*J12*J33+J21*J13*J32+J31*J12*J23-J31*J13*J22) *w2e-. .. 
(J12*J33-J13*J32)/(J11*J22*J33-J11*J23*J32-J21*J12*J33+. 
J21*J13*J32+J31*J12*J23-J31*J13*J22)*wle)*J32) *w2e+. .. 

(((J12*J33-J13*J32)/(J11*J22*J33-J11*J23*J32-. .. 
J21*J12*J33+J21*J13*J32+J31*J12*J23-J31*J13*J22)*w3e+... 
(J12*J23-J13*J22)/(J11*J22*J33-J11*J23*J32-J21*J12*J33+ . 
J21*J13*J32+J31*J12*J23-J31*J13*J22)*w2e) *J13+. .. 

((J22*J33-J23*J32)/(Jll*J22*J33-J11*J23*J32-J21*J12*J33+ 
J21*J13*J32+J31*J12*J23-J31*J13*J22) *w3e-. .. 

(J12*J23-J13*J22)/(J11*J22*J33-J11*J23*J32-J21*J12*J33+. 
J21*J13*J32+J31*J12*J23-J31*J13*J22)*wle) *J23+. .. 

(-(J22*J33-J23*J32)/(Jll*J22*J33-J11*J23*J32- ... 
J21*J12*J33+J21*J13*J32+J31*J12*J23-J31*J13*J22) *w2e-. .. 
(J12*J33-J13*J32)/(J11*J22*J33-J11*J23*J32-J21*J12*J33+ . 
J21*J13*J32+J31*J12*J23-J31*J13*J22)*wle)*J33)*w3e)+... 
(-1/2*al*q3e+((J12*J23-J13*J22)/(Jll*J22*J33- ... 
J11*J23*J32-J21*J12*J33+J21*J13*J32+J31*J12*J23-. .. 

J31*J13*J22)*J11-(J22*J33-J23*J32)/(Jll*J22*J33- ... 
J11*J23*J32-J21*J12*J33+J21*J13*J32+J31*J12*J23-. .. 

J31*J13*J22)*J31)*wle+J23-J13*J22)/(Jll*J22*J33- . . . 
J11*J23*J32-J21*J12*J33+J21*J13*J32+J31*J12*J23-. .. 

J31*J13*J22)*J12-(J22*J33-J23*J32)/ ... 

(J11*J22*J33-J11*J23*J32-J21*J12*J33+J21*J13*J32+... 
J31*J12*J23-J31*J13*J22)*J32)*w2e+((J12*J33-J13*J32)/ ... 
(J11*J22*J33-J11*J23*J32-J21*J12*J33+J21*J13*J32+. . . 
J31*J12*J23-J31*J13*J22)*w3e+(J12*J23-J13*J22) / . . . 

(J11*J22*J33-J11*J23*J32-J21*J12*J33+J21*J13*J32+... 
J31*J12*J23-J31*J13*J22)*w2e)*J12+((J22*J33-J23*J32) / . . . 
(J11*J22*J33-J11*J23*J32-J21*J12*J33+J21*J13*J32+... 
J31*J12*J23-J31*J13*J22)*w3e-(J12*J23-J13*J22)/ ... 

(J11*J22*J33-J11*J23*J32-J21*J12*J33+J21*J13*J32+. . . 
J31*J12*J23-J31*J13*J22)*wle)*J22+(-(J22*J33-J23*J32)/ .. 
(J11*J22*J33-J11*J23*J32-J21*J12*J33+J21*J13*J32+. . . 
J31*J12*J23-J31*J13*J22)*w2e-(J12*J33-J13*J32)/ ... 

(J11*J22*J33-J11*J23*J32-J21*J12*J33+J21*J13*J32+... 
J31*J12*J23-J31*J13*J22)*wle)*J32+((J12*J23-J13*J22) / . . . 
(J11*J22*J33-J11*J23*J32-J21*J12*J33+J21*J13*J32+... 
J31*J12*J23-J31*J13*J22)*J13-(J22*J33-J23*J32)/ ... 

(J11*J22*J33-J11*J23*J32-J21*J12*J33+J21*J13*J32+. . . 
J31*J12*J23-J31*J13*J22)*J33)*w3e)*(((-(Jll*J33- ... 

J13*J31)/(J11*J22*J33-J11*J23*J32-J21*J12*J33+. .. 
J21*J13*J32+J31*J12*J23-J31*J13*J22)*w3e- (J11*J23-. .. 
J13*J21)/(J11*J22*J33-J11*J23*J32-J21*J12*J33+. .. 
J21*J13*J32+J31*J12*J23-J31*J13*J22)*w2e) *J11+. .. 

(-(J21*J33-J23*J31)/(J11*J22*J33-J11*J23*J32-. .. 
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J21*J12*J33+J21*J13*J32+J31*J12*J23-J31*J13*J22) *w3e+. .. 

(J11*J23-J13*J21)/(J11*J22*J33-J11*J23*J32-J21*J12*J33+ . . . 
J21*J13*J32+J31*J12*J23-J31*J13*J22)*wle) *J21+. .. 

((J21*J33-J23*J31)/(J11*J22*J33-J11*J23*J32-J21*J12*J33+ .. 
J21*J13*J32+J31*J12*J23-J31*J13*J22)*w2e+ (J11*J33-. . . 
J13*J31)/(J11*J22*J33-J11*J23*J32-J21*J12*J33+ . . . 
J21*J13*J32+J31*J12*J23-J31*J13*J22)*wle)*J31) *wle+. .. 

((-(J11*J33-J13*J31)/(J11*J22*J33-J11*J23*J32~. . . 
J21*J12*J33+J21*J13*J32+J31*J12*J23-J31*J13*J22)*w3e- ... 

(J11*J23-J13*J21)/(J11*J22*J33-J11*J23*J32-J21*J12*J33+... 
J21*J13*J32+J31*J12*J2 3-J31*J13*J22)*w2e)* J12+. . . 

(-(J21*J33-J23*J31)/(J11*J22*J33-J11*J23*J32~. . . 
J21*J12*J33+J21*J13*J32+J31*J12*J23-J31*J13*J22) *w3e+. .. 

(J11*J23-J13*J21)/(J11*J22*J33-J11*J23*J32-J21*J12*J33+... 
J21*J13*J32+J31*J12*J23-J31*J13*J22)*wle)*J22+ ... 

((J21*J33-J23*J31)/(J11*J22*J33-J11*J23*J32-. .. 
J21*J12*J33+J21*J13*J32+J31*J12*J23-J31*J13*J22)*w2e+... 

(J11*J33-J13*J31)/(J11*J22*J33-J11*J23*J32-J21*J12*J33+... 
J21*J13*J32+J31*J12*J23-J31*J13*J22)*wle)*J32)*w2e+... 

((-(J11*J33-J13*J31)/(J11*J22*J33-J11*J23*J32-. . . 
J21*J12*J33+J21*J13*J32+J31*J12*J23-J31*J13*J22) *w3e-. .. 

(J11*J23-J13*J21)/(J11*J22*J33-J11*J23*J32-J21*J12*J33+ . . . 
J21*J13*J32+J31*J12*J23-J31*J13*J22)*w2e)* J13+. .. 

(-(J21*J33-J23*J31)/(J11*J22*J33-J11*J23*J32-. .. 
J21*J12*J33+J21*J13*J32+J31*J12*J23-J31*J13*J22)*w3e+. .. 

(J11*J23-J13*J21)/(J11*J22*J33-J11*J23*J32-J21*J12*J33+. 
J21*J13*J32+J31*J12*J23-J31*J13*J22)*wle)*J23+ ... 

((J21*J33-J23*J31)/(J11*J22*J33-J11*J23*J32-J21*J12*J33+ 
J21*J13*J32+J31*J12*J23-J31*J13*J22)*w2e+. .. 

(J11*J33-J13*J31)/(J11*J22*J33-J11*J23*J32-J21*J12*J33+. 
J21*J13*J32+J31*J12*J23-J31*J13*J22)*wle)*J33)*w3e)+... 
(l/2*al*q2e+((J12*J33-J13*J32)/(J11*J22*J33-J11*J23*J32- 
J21*J12*J33+J21*J13*J32+J31*J12*J23-J31*J13*J22)*J11+. . . 

(J22*J33-J23*J32)/(Jll*J22*J33-J11*J23*J32-J21*J12*J33+. 
J21*J13*J32+J31*J12*J23-J31*J13*J22)*J21)*wle+ ... 

((J12*J33-J13*J32)/(J11*J22*J33-J11*J23*J32-J21*J12*J33+ 
J21*J13*J32+J31*J12*J23-J31*J13*J22)*J12+ ... 

(J22*J33-J23*J32)/(Jll*J22*J33-J11*J23*J32-J21*J12*J33+. 
J21*J13*J32+J31*J12*J23-J31*J13*J22)*J22)*w2e+ ... 

((J12*J33-J13*J32)/(J11*J22*J33-J11*J23*J32-J21*J12*J33+ 
J21*J13*J32+J31*J12*J23-J31*J13*J22)*J13+ . . . 

(J22*J33-J23*J32)/(Jll*J22*J33-J11*J23*J32-J21*J12*J33+. 
J21*J13*J32+J31*J12*J23-J31*J13*J22)*J23)*w3e+ ... 

((J12*J33-J13*J32)/(J11*J22*J33-J11*J23*J32-J21*J12*J33+ 
J21*J13*J32+J31*J12*J23-J31*J13*J22)*w3e+ ... 

(J12*J23-J13*J22)/(J11*J22*J33-J11*J23*J32-J21*J12*J33+. 
J21*J13*J32+J31*J12*J23-J31*J13*J22)*w2e)* J13+. .. 

((J22*J33-J23*J32)/(Jll*J22*J33-J11*J23*J32-J21*J12*J33+ 
J21*J13*J32+J31*J12*J23-J31*J13*J22) *w3e-. .. 

(J12*J23-J13*J22)/(J11*J22*J33-J11*J23*J32-J21*J12*J33+. 
J21*J13*J32+J31*J12*J23-J31*J13*J22)*wle)*J23+ ... 

(-(J22*J33-J23*J32)/(Jll*J22*J33-J11*J23*J32- ... 
J21*J12*J33+J21*J13*J32+J31*J12*J23-J31*J13*J22) *w2e-. .. 

(J12*J33-J13*J32)/(J11*J22*J33-J11*J23*J32-J21*J12*J33+... 
J21*J13*J32+J31*J12*J23-J31*J13*J22)*wle)*J33)* ... 
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((((J11*J32-J12*J31)/(J11*J22*J33-J11*J23*J32-. .. 
J21*J12*J33+J21*J13*J32+J31*J12*J23-J31*J13*J22) *w3e+. .. 
(J11*J22-J12*J21)/(J11*J22*J33-J11*J23*J32-J21*J12*J33+. 
J21*J13*J32+J31*J12*J23-J31*J13*J22)*w2e)*J11+((J21*J32- 
J22*J31)/(J11*J22*J33-J11*J23*J32-J21*J12*J33+ ... 
J21*J13*J32+J31*J12*J23-J31*J13*J22)*w3e- (J11*J22-. .. 
J12*J21)/(J11*J22*J33-J11*J23*J32-J21*J12*J33+ ... 
J21*J13*J32+J31*J12*J23-J31*J13*J22)*wle) *J21+. .. 

(-(J21*J32-J22*J31)/(J11*J22*J33-J11*J23*J32~. . . 
J21*J12*J33+J21*J13*J32+J31*J12*J23-J31*J13*J22) *w2e-. .. 
(J11*J32-J12*J31)/(J11*J22*J33-J11*J23*J32-J21*J12*J33+. 
J21*J13*J32+J31*J12*J23~. .. 

J31*J13*J22)*wle)*J31)*wle+(((Jll*J32-J12*J31)/ ... 

(J11*J22*J33-J11*J23*J32-J21*J12*J33+J21*J13*J32+. . . 
J31*J12*J23-J31*J13*J22)*w3e+(Jll*J22-J12*J21) / . . . 

(J11*J22*J33-J11*J23*J32-J21*J12*J33+J21*J13*J32+... 
J31*J12*J23-J31*J13*J22)*w2e)*J12+((J21*J32-J22*J31)/ . . . 
(J11*J22*J33-J11*J23*J32-J21*J12*J33+J21*J13*J32+ . . . 
J31*J12*J23-J31*J13*J22) *w3e- (Jll*J22-J12*J21)/ ... 

(J11*J22*J33-J11*J23*J32-J21*J12*J33+J21*J13*J32+... 
J31*J12*J2 3-J31*J13*J22)*wle)*J22+(-(J21*J32-J22*J31) / . . 
(J11*J22*J33-J11*J23*J32-J21*J12*J33+J21*J13*J32+... 
J31*J12*J23-J31*J13*J22)*w2e-(Jll*J32-J12*J31)/ ... 

(J11*J22*J33-J11*J23*J32-J21*J12*J33+J21*J13*J32+ . . . 
J31*J12*J23-J31*J13*J22)*wle)*J32)*w2e+(( (J11*J32-. .. 
J12*J31)/(J11*J22*J33-J11*J23*J32-J21*J12*J33+ ... 
J21*J13*J32+J31*J12*J23-J31*J13*J22)*w3e+ (J11*J22-. . . 
J12*J21)/(J11*J22*J33-J11*J23*J32-J21*J12*J33+ . . . 
J21*J13*J32+J31*J12*J23-J31*J13*J22)*w2e)* J13+. .. 

((J21*J32-J22*J31)/(J11*J22*J33-J11*J23*J32-J21*J12*J33+ 
J21*J13*J32+J31*J12*J23-J31*J13*J22)*w3e- (J11*J22-. .. 
J12*J21)/(J11*J22*J33-J11*J23*J32-J21*J12*J33+. . . 
J21*J13*J32+J31*J12*J23-J31*J13*J22)*wle) *J23+. .. 

(-(J21*J32-J22*J31)/(J11*J22*J33-J11*J23*J32-. .. 
J21*J12*J33+J21*J13*J32+J31*J12*J23-J31*J13*J22) *w2e-. .. 
(J11*J32-J12*J31)/(J11*J22*J33-J11*J23*J32-J21*J12*J33+. 
J21*J13*J32+J31*J12*J23-J31*J13*J22)*wle)*J33)*w3e); 
Lfphi = al*(l/2*q0e*wle-l/2*q3e*w2e+l/2*q2e*w3e)+(( (J12*J33-. .. 

J13*J32)/(J11*J22*J33-J11*J23*J32-J21*J12*J33+ ... 
J21*J13*J32+J31*J12*J23-J31*J13*J22)*w3e+. .. 

(J12*J23-J13*J22)/(J11*J22*J33-J11*J23*J32-. .. 
J21*J12*J33+J21*J13*J32+J31*J12*J23-. .. 

J31*J13*J22)*w2e)*J11+((J22*J33-J23*J32)/ ... 

(J11*J22*J33-J11*J23*J32-J21*J12*J33+J21*J13*J32+... 
J31*J12*J23-J31*J13*J22)*w3e-(J12*J2 3-J13*J22) / . . . 

(J11*J22*J33-J11*J23*J32-J21*J12*J33+J21*J13*J32+ . . . 
J31*J12*J23-J31*J13*J22)*wle)*J21+(-(J22*J33-J23*J32)/ .. 
(J11*J22*J33-J11*J23*J32-J21*J12*J33+J21*J13*J32+... 
J31*J12*J23-J31*J13*J22) *w2e- (J12*J33-J13*J32)/ ... 

(J11*J22*J33-J11*J23*J32-J21*J12*J33+J21*J13*J32+... 
J31*J12*J2 3-J31*J13*J22)*wle)*J31)*wle+(( (J12*J33-. . . 
J13*J32)/(J11*J22*J33-J11*J23*J32-J21*J12*J33+. .. 
J21*J13*J32+J31*J12*J23-J31*J13*J22) *w3e+. .. 

(J12*J23-J13*J22)/(J11*J22*J33-J11*J23*J32-J21*J12*J33+. 


151 


J21*J13*J32+J31*J12*J23-J31*J13*J22)*w2e) *J12+. .. 

((J22*J33-J23*J32)/(Jll*J22*J33-J11*J2 3*J32-J21*J12*J33+ . . . 
J21*J13*J32+J31*J12*J23-J31*J13*J22)*w3e- (J12*J23~. .. 

J13*J22)/(J11*J22*J33-J11*J23*J32-J21*J12*J33+. .. 
J21*J13*J32+J31*J12*J23-J31*J13*J22)*wle) *J22+. .. 

(-(J22*J33-J23*J32)/(Jll*J22*J33-J11*J23*J32- . . . 
J21*J12*J33+J21*J13*J32+J31*J12*J23-J31*J13*J22) *w2e-. .. 

(J12*J33-J13*J32)/(J11*J22*J33-J11*J23*J32-. . . 
J21*J12*J33+J21*J13*J32+J31*J12*J23-. .. 

J31*J13*J22)*wle)*J32)*w2e+(((J12*J33-J13*J32) / . . . 
(J11*J22*J33-J11*J2 3*J32-J21*J12*J33+J21*J13*J32+. . . 
J31*J12*J23-J31*J13*J22)*w3e+(J12*J23-J13*J22) / . . . 

(J11*J22*J33-J11*J23*J32-J21*J12*J33+J21*J13*J32+ ... 
J31*J12*J23-J31*J13*J22)*w2e)*J13+( (J22*J33-. .. 

J23*J32)/(J11*J22*J33-J11*J23*J32-J21*J12*J33+. .. 
J21*J13*J32+J31*J12*J23-J31*J13*J22)*w3e- (J12*J23-. . . 

J13*J22)/(J11*J22*J33-J11*J23*J32-J21*J12*J33+. .. 
J21*J13*J32+J31*J12*J23-J31*J13*J22)*wle) *J23+. .. 

(-(J22*J33-J23*J32)/(Jll*J22*J33-J11*J23*J32- ... 
J21*J12*J33+J21*J13*J32+J31*J12*J23-. .. 

J31*J13*J22)*w2e-(J12*J33-J13*J32)/(Jll*J22*J33- ... 
J11*J23*J32-J21*J12*J33+J21*J13*J32+. .. 
J31*J12*J23-J31*J13*J22)*wle)*J33)*w3e; 


phidot = 
b 

wu = 

wa = 

qa = 

%qu 


wle + al*qle; 

Lfphi; 

-Lf2phi - 2*gamma*Lfphi - gamma A 2*phi; 
wle; 

[w2e;w3e]; 
qve(2:3); 


Fw 

= 

-inv(J) 

*wx 

S21 

= 

Fw(2:3, 

l); 

S22 


Fw(2:3, 

2:3 

if 

' (a 

'*a) >= 

Bet 


adinv 

= a/ ( 

a' * 


else adinv = a/Betal A 2; 
end 

Pd = eye(2) - adinv*a'; 

udbar = adinv*b; 
if select == 1 

%% Lyapunov Function-Based Control Law 
c = wa'*Pd*wa; 
if sqrt(c) >= Beta2 

etad = wa/(wa'*Pd*wa); 
else etad = wa/Beta2; 
end 

y = etad*(-wa'*S21*wu - wa'*S22*wa - wa'*udbar - k*qve'*we) - d*wa; 

else 

%% Feedback Linearizing Control Law 
y = (-S21*wu - S22*wa - k*qa - d*wa); 


end 
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tau_a 

tau 

T 


= udbar + Pd*y 
= [0;tau_a]; 

= J*tau; 
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